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

[Romaric Jodin via cfe-commits]

https://github.com/rjodinchr created https://github.com/llvm/llvm-project/pull/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.

>From aed077ba203b15df301e81a637641ea74c7f138c Mon Sep 17 00:00:00 2001
From: Romaric Jodin <rjodin at google.com>
Date: Tue, 20 Oct 2026 11:13:49 +0200
Subject: [PATCH] [libclc] Support fp32 subnormal inputs and outputs in SPIR-V
 sw fma

`__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.

Change-Id: I7368ef462b16bdc5f0c05704052cb9e2feb86ea2
---
 libclc/clc/lib/spirv/math/clc_sw_fma.cl | 50 +++++++++++++++++--------
 1 file changed, 35 insertions(+), 15 deletions(-)

diff --git a/libclc/clc/lib/spirv/math/clc_sw_fma.cl b/libclc/clc/lib/spirv/math/clc_sw_fma.cl
index e73f53e6bcd769..fae76e62050eb8 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));
 }