[AGC] Emit v_pk_fma_f16 with exact single rounding (#420)

Completes the fused-FMA slice deferred by the VOP3P first slice (#145).
v_pk_fma_f16 previously failed emission loudly because an f32
multiply-add followed by an f16 pack rounds twice; the pinned miss is
fma(0x4100, 0x7522, 0x04EA) = 0x7A6B fused vs 0x7A6A via f32.

The f32 product of two f16 values is exact, so only the addition needs
correcting: compute sum = RN(product + addend), recover the exact
residual with Knuth 2Sum, and if the sum is inexact with an even
significand, step one ulp towards the true value. That is round-to-odd,
and rounding the f32 result to f16 with round-to-nearest-even then
matches a true fused f16 FMA exactly (24 significand bits >= 11 + 2).
Inf/NaN inputs turn the residual into NaN, the ordered compare skips the
parity fix, and IEEE special behaviour passes through unchanged. The
op_sel/op_sel_hi/neg_lo/neg_hi source modifiers apply to src2 through
the existing operand path; clamp stays rejected like the other packed
ops.

Every op in the 2Sum chain is decorated NoContraction: without it the
AMD RDNA3 Windows driver folds the sequence, collapses the residual to
zero, and the midpoint case decays to the double-rounded result. This
was caught by running the emitted shader on a real device (see below).

Verification:
- A mirror of the emitted sequence was checked against an exact
  integer reference (every finite f16 is m * 2^-24, so a*b + c is an
  exact Int128 multiple of 2^-48, rounded once to f16 RNE) across 34M
  cases: directed midpoint pins, random sweeps over all operand
  classes, tiny-addend midpoint stress, subnormal products, and
  Inf/NaN propagation. 0 mismatches.
- ShaderDump gains a pk-f16 program covering all five packed opcodes,
  both fma modifier paths, and the pinned constants; all programs
  decode and emit.
- The executable exec program now computes the pinned fma and its
  negated-addend twin (0x7A6B7A6B / 0x7A6A7A6A, straddling an f16
  midpoint) and stores them at offsets 20/24; GpuConformance checks
  both on device. All values match on an AMD Radeon RX 7700 XT.

src/SharpEmu.ShaderCompiler.Vulkan/Gen5SpirvTranslator.Alu.cs

@@ -953,22 +953,13 @@ public static partial class Gen5SpirvTranslator
                 case "VPkMulF16":
                 case "VPkMinF16":
                 case "VPkMaxF16":
+                case "VPkFmaF16":
                     if (!TryEmitPackedF16(instruction, out result, out error))
                     {
                         return false;
                     }
 
                     break;
-                case "VPkFmaF16":
-                    // Deliberately loud: a fused f16 FMA rounds the product+add once,
-                    // whereas doing the multiply-add in f32 and rounding to f16 at the
-                    // end double-rounds. Concrete miss: fma(0x4100, 0x7522, 0x04EA) is
-                    // 0x7A6B fused but 0x7A6A via f32. Exact emulation (round-to-odd
-                    // f32 product then RNE pack) is a planned follow-up slice.
-                    error =
-                        $"unsupported vop3p opcode {instruction.Opcode} " +
-                        "(fused f16 FMA requires single-rounding; deferred to a later slice)";
-                    return false;
                 default:
                     error = $"unsupported vector opcode {instruction.Opcode}";
                     return false;
@@ -1008,8 +999,9 @@ public static partial class Gen5SpirvTranslator
         // even. For add and mul this is bit-exact to a true f16 op (the f32 result
         // rounds losslessly to f16 by the double-rounding theorem; a f16 product even
         // fits in f32 exactly). min/max carry no rounding, so they are exact once the
-        // conversions are. v_pk_fma_f16 is intentionally not routed here because a
+        // conversions are. v_pk_fma_f16 cannot be reproduced by a plain f32
-        // fused f16 FMA cannot be reproduced by an f32 multiply-add plus a pack.
+        // multiply-add plus a pack (that double-rounds), so it goes through the
+        // round-to-odd sequence in EmitPackedF16FusedMultiplyAdd instead.
         private bool TryEmitPackedF16(
             Gen5ShaderInstruction instruction,
             out uint result,
@@ -1029,7 +1021,8 @@ public static partial class Gen5SpirvTranslator
                 return false;
             }
 
-            for (var index = 0; index < 2; index++)
+            var sourceCount = instruction.Opcode == "VPkFmaF16" ? 3 : 2;
+            for (var index = 0; index < sourceCount; index++)
             {
                 var source = instruction.Sources[index];
                 if (source.Kind is not (Gen5OperandKind.VectorRegister or Gen5OperandKind.ScalarRegister))
@@ -1054,6 +1047,12 @@ public static partial class Gen5SpirvTranslator
         {
             var left = EmitPackedF16Operand(instruction, control, 0, highLane);
             var right = EmitPackedF16Operand(instruction, control, 1, highLane);
+            if (instruction.Opcode == "VPkFmaF16")
+            {
+                var addend = EmitPackedF16Operand(instruction, control, 2, highLane);
+                return EmitFloatToHalf(EmitPackedF16FusedMultiplyAdd(left, right, addend));
+            }
+
             var value = instruction.Opcode switch
             {
                 "VPkAddF16" => _module.AddInstruction(SpirvOp.FAdd, _floatType, left, right),
@@ -1065,6 +1064,75 @@ public static partial class Gen5SpirvTranslator
             return EmitFloatToHalf(Bitcast(_uintType, value));
         }
 
+        // Fused f16 multiply-add with a single rounding, emulated in f32 without the
+        // Float16 capability. The f32 product of two widened f16 values is exact
+        // (11-bit significands, and the exponent stays inside the f32 normal range:
+        // any non-zero product magnitude is in [2^-48, 2^33]), so only the addition
+        // rounds. An f32 add then an f16 pack would round twice; instead the add is
+        // corrected to round-to-odd, which a following round-to-nearest-even pack
+        // turns into the exactly-once-rounded fused result (innocuous double rounding
+        // holds because f32 carries 24 significand bits >= 11 + 2).
+        //
+        // sum = RN(product + addend); Knuth's 2Sum recovers the exact residual
+        // (product + addend) - sum from four more RN ops. 2Sum is exact for any two
+        // finite f32 inputs; no intermediate here can overflow (|product| < 2^33,
+        // |addend| < 2^16) and none can enter the f32 subnormal range (every finite
+        // value in play is a multiple of 2^-48 by construction), so implementation
+        // f32 denorm-flush modes never see a denormal. If the residual says the sum
+        // was inexact and the sum's significand is even, step one ulp towards the
+        // true value: consecutive floats have consecutive sign-magnitude encodings,
+        // so that neighbour is the enclosing float with the odd significand.
+        //
+        // Inf/NaN inputs make the residual NaN (e.g. sum - addend = Inf - Inf); the
+        // ordered compare below is then false and the IEEE sum passes through
+        // unchanged. A residual of zero also covers the exact-sum case, where the
+        // parity fix must not fire. Returns the round-to-odd f32 bit pattern.
+        private uint EmitPackedF16FusedMultiplyAdd(uint left, uint right, uint addend)
+        {
+            var product = EmitPreciseFloat(SpirvOp.FMul, left, right);
+            var sum = EmitPreciseFloat(SpirvOp.FAdd, product, addend);
+
+            var productPart = EmitPreciseFloat(SpirvOp.FSub, sum, addend);
+            var addendPart = EmitPreciseFloat(SpirvOp.FSub, sum, productPart);
+            var productError = EmitPreciseFloat(SpirvOp.FSub, product, productPart);
+            var addendError = EmitPreciseFloat(SpirvOp.FSub, addend, addendPart);
+            var residual = EmitPreciseFloat(SpirvOp.FAdd, productError, addendError);
+
+            var sumBits = Bitcast(_uintType, sum);
+            var residualBits = Bitcast(_uintType, residual);
+            var inexact = _module.AddInstruction(
+                SpirvOp.FOrdNotEqual, _boolType, residual, Float(0));
+            var evenSignificand = Equal(BitwiseAnd(sumBits, UInt(1)), 0);
+            var adjust = _module.AddInstruction(
+                SpirvOp.LogicalAnd, _boolType, inexact, evenSignificand);
+
+            // Residual sign relative to the sum picks the step direction: same sign
+            // means the true value lies away from zero (encoding + 1), opposite sign
+            // means towards zero (encoding - 1). The sum cannot be zero here (any
+            // inexact sum has magnitude >= 2^-48) and cannot be the largest finite
+            // value (its significand is odd), so the step never crosses zero or Inf.
+            var towardZero = IsNotZero(
+                BitwiseAnd(BitwiseXor(sumBits, residualBits), UInt(0x8000_0000)));
+            var stepped = SelectU(
+                towardZero,
+                ISubU(sumBits, UInt(1)),
+                IAdd(sumBits, UInt(1)));
+            return SelectU(adjust, stepped, sumBits);
+        }
+
+        // A float op the driver must evaluate exactly as written. The 2Sum
+        // residual above is error-free only op by op; without NoContraction
+        // driver compilers fold the sequence (e.g. contract product+sum into an
+        // f32 fma and simplify the rebuilt terms), collapsing the residual to
+        // zero. Observed on AMD RDNA3 Windows: the pinned midpoint case decays
+        // to the double-rounded result unless every op in the chain is marked.
+        private uint EmitPreciseFloat(SpirvOp operation, uint left, uint right)
+        {
+            var value = _module.AddInstruction(operation, _floatType, left, right);
+            _module.AddDecoration(value, SpirvDecoration.NoContraction);
+            return value;
+        }
+
         // Reads source `index`, selects the half feeding this lane (op_sel / op_sel_hi),
         // widens it exactly to f32 and applies the lane's negate modifier (neg_lo / neg_hi).
         private uint EmitPackedF16Operand(

src/SharpEmu.ShaderCompiler.Vulkan/SpirvModuleBuilder.cs

@@ -238,6 +238,7 @@ public enum SpirvDecoration : uint
     Binding = 33,
     DescriptorSet = 34,
     Offset = 35,
+    NoContraction = 42,
 }
 
 public enum SpirvBuiltIn : uint

tools/SharpEmu.Tools.GpuConformance/Program.cs

@@ -12,6 +12,11 @@
 //   [2] v_mul_lo_i32 -> low  32 bits of the same product
 //   [3] store attempted with EXEC=0 -> must NOT land (sentinel remains)
 //   [4] store after EXEC restored   -> 1.5f (0x3FC00000)
+//   [5] v_pk_fma_f16 fma(2.5h, 21024h,  7.496e-5h) -> 0x7A6B packed; the exact
+//       sum sits just above an f16 midpoint, so a double-rounded f32
+//       multiply-add would give 0x7A6A instead
+//   [6] the same fma with the addend negated -> 0x7A6A packed (just below the
+//       same midpoint), pinning the opposite rounding direction
 // Every other word of the buffer must still hold the sentinel afterwards.
 //
 // Creating the compute pipeline doubles as a driver-acceptance check for the
@@ -35,6 +40,12 @@ var expectedHi = (uint)(product >> 32);
 var expectedLo = (uint)product;
 var expectedRestored = BitConverter.SingleToUInt32Bits(1.5f);
 
+// v_pk_fma_f16 of (0x4100, 0x7522, 0x04EA) per lane: the exact product
+// 2.5 * 21024 = 52560 is an f16 tie (between 0x7A6A and 0x7A6B), so the tiny
+// addend decides the rounding direction under a single fused rounding.
+const uint ExpectedPkFma = 0x7A6B_7A6B;
+const uint ExpectedPkFmaNeg = 0x7A6A_7A6A;
+
 unsafe
 {
     var spvPath = args.Length > 0
@@ -352,6 +363,8 @@ unsafe
         ("v_mul_lo_i32 lo(0x7FFFFFFF*0x10003)", words[2], expectedLo),
         ("exec=0 store suppressed (offset 12 sentinel)", words[3], Sentinel),
         ("store after exec restore (offset 16)", words[4], expectedRestored),
+        ("v_pk_fma_f16 fused rounds up at midpoint", words[5], ExpectedPkFma),
+        ("v_pk_fma_f16 neg addend rounds down", words[6], ExpectedPkFmaNeg),
     };
     var failures = 0;
     foreach (var (name, actual, expected) in results)

tools/SharpEmu.Tools.ShaderDump/Program.cs

@@ -42,6 +42,23 @@ const ulong ProgramAddress = 0x100000;
         0xD56C0005, 0x00020501, // v_mul_hi_i32 v5, v1, v2
         0xBF810000,             // s_endpgm
     ]),
+    // Packed f16 (VOP3P) arithmetic, including the fused multiply-add. The
+    // constants pin the double-rounding regression from the VOP3P first slice:
+    // fma(0x4100, 0x7522, 0x04EA) must round once to 0x7A6B (an f32
+    // multiply-add then pack yields 0x7A6A). The last fma exercises the src2
+    // neg_lo/neg_hi modifier path.
+    ("pk-f16", true, [
+        0x7E0002FF, 0x41004100, // v_mov_b32 v0, 0x41004100 (2.5 packed)
+        0x7E0202FF, 0x75227522, // v_mov_b32 v1, 0x75227522 (21024 packed)
+        0x7E0402FF, 0x04EA04EA, // v_mov_b32 v2, 0x04EA04EA (~7.496e-5 packed)
+        0xCC0E4003, 0x1C0A0300, // v_pk_fma_f16 v3, v0, v1, v2
+        0xCC0F4004, 0x18020500, // v_pk_add_f16 v4, v0, v2
+        0xCC104005, 0x18020300, // v_pk_mul_f16 v5, v0, v1
+        0xCC114006, 0x18020300, // v_pk_min_f16 v6, v0, v1
+        0xCC124007, 0x18020300, // v_pk_max_f16 v7, v0, v1
+        0xCC0E4408, 0x9C0A0300, // v_pk_fma_f16 v8, v0, v1, neg_lo:[0,0,1] neg_hi:[0,0,1] v2
+        0xBF810000,             // s_endpgm
+    ]),
     ("mrt", true, [
         0x7E0002FF, 0x3F800000, // v_mov_b32 v0, 1.0f
         0x7E0202FF, 0x00000000, // v_mov_b32 v1, 0.0f
@@ -115,7 +132,10 @@ const ulong ProgramAddress = 0x100000;
     // Executable end-to-end test: compute with real ALU instructions, then
     // buffer_store_dword results to guestBuffers[0] at offsets 0/4/8, prove
     // that a store with EXEC=0 does not land (offset 12 stays sentinel), and
-    // that stores work again after EXEC is restored (offset 16).
+    // that stores work again after EXEC is restored (offset 16). Offsets 20/24
+    // hold the packed fused f16 FMA and its negated-addend twin, whose exact
+    // results (0x7A6B7A6B / 0x7A6A7A6A) straddle an f16 midpoint and therefore
+    // catch any double-rounding regression on real hardware.
     ("exec", true, [
         0xBFA10001,             // s_clause 0x1 (hint no-op in an executed program, needs #108)
         0x7E0002FF, 0x3FC00000, // v_mov_b32 v0, 1.5f
@@ -133,6 +153,13 @@ const ulong ProgramAddress = 0x100000;
         0xE070000C, 0x80020200, // buffer_store_dword v2, off, s[8:11], 0 offset:12 (masked, must not land)
         0xBEFE03C1,             // s_mov_b32 exec_lo, -1      -> lane active again
         0xE0700010, 0x80020000, // buffer_store_dword v0, off, s[8:11], 0 offset:16
+        0x7E0E02FF, 0x41004100, // v_mov_b32 v7, 0x41004100 (2.5 packed)
+        0x7E1002FF, 0x75227522, // v_mov_b32 v8, 0x75227522 (21024 packed)
+        0x7E1202FF, 0x04EA04EA, // v_mov_b32 v9, 0x04EA04EA (~7.496e-5 packed)
+        0xCC0E400A, 0x1C261107, // v_pk_fma_f16 v10, v7, v8, v9
+        0xCC0E440B, 0x9C261107, // v_pk_fma_f16 v11, v7, v8, neg_lo:[0,0,1] neg_hi:[0,0,1] v9
+        0xE0700014, 0x80020A00, // buffer_store_dword v10, off, s[8:11], 0 offset:20
+        0xE0700018, 0x80020B00, // buffer_store_dword v11, off, s[8:11], 0 offset:24
         0xBF810000,             // s_endpgm
     ]),
 ];