[libclc] c6600a4 - [libclc] Support fp32 subnormal inputs and outputs in SPIR-V sw fma (#233791)

[via cfe-commits]

Author: Romaric Jodin
Date: 2026-10-26T06:42:21+08:00
New Revision: c6600a4a5b4037af9fe38e2d93144c2c1f1f2818

URL: https://github.com/llvm/llvm-project/commit/c6600a4a5b4037af9fe38e2d93144c2c1f1f2818
DIFF: https://github.com/llvm/llvm-project/commit/c6600a4a5b4037af9fe38e2d93144c2c1f1f2818.diff

LOG: [libclc] Support fp32 subnormal inputs and outputs in SPIR-V sw fma (#233791)

`__clc_sw_fma` previously assumed all non-zero inputs were normal floats
(extracting biased exponent `0` as `-127` with the implicit leading bit
`0x800000` set) and unconditionally flushed results with exponent `<=
-127` to zero. On devices that support fp32 subnormals and use software
fma, this caused `fma`, `expm1`, and `erf` to fail the CTS.

Normalize non-zero subnormal inputs using `__clc_clz` and, when fp32
subnormals are supported (`!__clc_denormals_are_zero_fp32()`),
right-shift the significand before rounding for exponents in `[-150,
-127]` so subnormal outputs and subnormal-to-normal rounding overflows
are rounded to nearest-even.

Added: 
    

Modified: 
    libclc/clc/lib/spirv/math/clc_sw_fma.cl

Removed: 
    


################################################################################
diff  --git a/libclc/clc/lib/spirv/math/clc_sw_fma.cl b/libclc/clc/lib/spirv/math/clc_sw_fma.cl
index e73f53e6bcd769..c2f67a80e2a668 100644
--- a/libclc/clc/lib/spirv/math/clc_sw_fma.cl
+++ b/libclc/clc/lib/spirv/math/clc_sw_fma.cl
@@ -142,16 +142,24 @@ _CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) {
 
   struct fp st_a, st_b, st_c;
 
-  st_a.exponent = a == .0f ? 0 : ((__clc_as_uint(a) & 0x7f800000) >> 23) - 127;
-  st_b.exponent = b == .0f ? 0 : ((__clc_as_uint(b) & 0x7f800000) >> 23) - 127;
-  st_c.exponent = c == .0f ? 0 : ((__clc_as_uint(c) & 0x7f800000) >> 23) - 127;
+  uint ua = __clc_as_uint(a) & 0x7fffffff;
+  uint ub = __clc_as_uint(b) & 0x7fffffff;
+  uint uc = __clc_as_uint(c) & 0x7fffffff;
+
+  uint sa = ua < 0x800000u ? __clc_clz(ua) - 8 : 0;
+  uint sb = ub < 0x800000u ? __clc_clz(ub) - 8 : 0;
+  uint sc = uc < 0x800000u ? __clc_clz(uc) - 8 : 0;
+
+  st_a.exponent = ua < 0x800000u ? -126 - (int)sa : (int)(ua >> 23) - 127;
+  st_b.exponent = ub < 0x800000u ? -126 - (int)sb : (int)(ub >> 23) - 127;
+  st_c.exponent = uc < 0x800000u ? -126 - (int)sc : (int)(uc >> 23) - 127;
 
   st_a.mantissa =
-      u2_set_u(a == .0f ? 0 : (__clc_as_uint(a) & 0x7fffff) | 0x800000);
+      u2_set_u(ua < 0x800000u ? (ua << sa) : (ua & 0x7fffff) | 0x800000);
   st_b.mantissa =
-      u2_set_u(b == .0f ? 0 : (__clc_as_uint(b) & 0x7fffff) | 0x800000);
+      u2_set_u(ub < 0x800000u ? (ub << sb) : (ub & 0x7fffff) | 0x800000);
   st_c.mantissa =
-      u2_set_u(c == .0f ? 0 : (__clc_as_uint(c) & 0x7fffff) | 0x800000);
+      u2_set_u(uc < 0x800000u ? (uc << sc) : (uc & 0x7fffff) | 0x800000);
 
   st_a.sign = __clc_as_uint(a) & 0x80000000;
   st_b.sign = __clc_as_uint(b) & 0x80000000;
@@ -214,6 +222,10 @@ _CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) {
     st_fma.sign = st_mul.sign ^ 0x80000000;
   }
 
+  if (u2_zero(st_fma.mantissa)) {
+    return 0.0f;
+  }
+
   // detect overflow/underflow
   int overflow_bits = 3 - u2_clz(st_fma.mantissa);
 
@@ -226,6 +238,19 @@ _CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) {
     overflow_bits = 0;
   }
 
+  if (st_fma.exponent <= -127) {
+    if (__clc_denormals_are_zero_fp32() || st_fma.exponent < -150) {
+      return __clc_as_float(st_fma.sign);
+    }
+    int sub_shift = -126 - st_fma.exponent;
+    uint2 sub_mask =
+        u2_add(u2_sll(u2_set_u(1), sub_shift), u2_set(0xffffffff, 0xffffffff));
+    cutoff_bits =
+        u2_or(cutoff_bits, !u2_zero(u2_and(st_fma.mantissa, sub_mask)));
+    st_fma.mantissa = u2_srl(st_fma.mantissa, sub_shift);
+    st_fma.exponent = -127;
+  }
+
   // rounding
   uint2 trunc_mask = u2_add(u2_sll(u2_set_u(1), C_ADJUST + overflow_bits),
                             u2_set(0xffffffff, 0xffffffff));
@@ -245,26 +270,21 @@ _CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) {
   // Shift mantissa back to bit 23
   st_fma.mantissa = u2_srl(st_fma.mantissa, C_ADJUST + overflow_bits);
 
+  if (st_fma.exponent == -127) {
+    return __clc_as_float(st_fma.sign | (uint)st_fma.mantissa.lo);
+  }
+
   // Detect rounding overflow
   if (u2_gt(st_fma.mantissa, u2_set_u(0xffffff))) {
     ++st_fma.exponent;
     st_fma.mantissa = u2_srl(st_fma.mantissa, 1);
   }
 
-  if (u2_zero(st_fma.mantissa)) {
-    return 0.0f;
-  }
-
   // Flating point range limit
   if (st_fma.exponent > 127) {
     return __clc_as_float(__clc_as_uint(INFINITY) | st_fma.sign);
   }
 
-  // Flush denormals
-  if (st_fma.exponent <= -127) {
-    return __clc_as_float(st_fma.sign);
-  }
-
   return __clc_as_float(st_fma.sign | ((st_fma.exponent + 127) << 23) |
                         ((uint)st_fma.mantissa.lo & 0x7fffff));
 }