[libclc] [libclc] Produce subnormal results for fmod on FP32 (PR #215768)

[via cfe-commits]

https://github.com/Lurie97 created https://github.com/llvm/llvm-project/pull/215768

The FP32 __clc_fmod reduction reconstructs the result's exponent with as_float(ey << 23) and rescales the reduced mantissa by it. When the divisor y is subnormal, ey == 0, so that factor is 0.0f and any result whose magnitude is subnormal (all results of fmod(x, subnormal_y)) is flushed to zero. On devices that report CL_FP_DENORM this loses the correct subnormal result (OpenCL-CTS math_brute_force fmod failed with the reference subnormal vs a returned 0, e.g.
fmod(-0x1.bcb23p-71, -0x1.8c9104p-127) -> 0 instead of -0x1.f79c38p-128).

Split the existing body into a normal-only core __clc_fmod_impl and wrap it: when subnormals are supported and y is subnormal, scale the operands up by 2^25 (turning the largest subnormal normal) so the core only sees a normal-sized modulus, then scale the result back down with __clc_ldexp, which produces subnormals naturally. Uses fmod(a*t, b*t) == t*fmod(a,b). For a normal x the huge value is first reduced modulo the scaled-up ys = y*2^25 (a multiple of y, so fmod(x, y) == fmod(x mod ys, y)) to avoid overflowing when scaling x up. A subnormal x with a normal y is unaffected: the core takes its ex < ey path and returns x directly.

On flush-to-zero devices __clc_denormals_are_zero_fp32() is true and the wrapper is skipped, preserving the previous behavior.

Verified against glibc with a scalar model over 40M random pairs with heavy subnormal coverage (x-sub, y-sub, both): zero mismatches. On i.MX95 Mali-G310 (Panfrost/Rusticl) advertising CL_FP_DENORM, every previously failing subnormal point now returns the correct result.

>From f3714b336a451af9270ad4f57a062ddb03401c6c Mon Sep 17 00:00:00 2001
From: jiajia Qian <jiajia.qian at nxp.com>
Date: Thu, 13 Aug 2026 16:10:03 +0800
Subject: [PATCH] [libclc] Produce subnormal results for fmod on FP32

The FP32 __clc_fmod reduction reconstructs the result's exponent with
as_float(ey << 23) and rescales the reduced mantissa by it. When the
divisor y is subnormal, ey == 0, so that factor is 0.0f and any result
whose magnitude is subnormal (all results of fmod(x, subnormal_y)) is
flushed to zero. On devices that report CL_FP_DENORM this loses the
correct subnormal result (OpenCL-CTS math_brute_force fmod failed with
the reference subnormal vs a returned 0, e.g.
fmod(-0x1.bcb23p-71, -0x1.8c9104p-127) -> 0 instead of -0x1.f79c38p-128).

Split the existing body into a normal-only core __clc_fmod_impl and wrap
it: when subnormals are supported and y is subnormal, scale the operands
up by 2^25 (turning the largest subnormal normal) so the core only sees
a normal-sized modulus, then scale the result back down with __clc_ldexp,
which produces subnormals naturally. Uses fmod(a*t, b*t) == t*fmod(a,b).
For a normal x the huge value is first reduced modulo the scaled-up
ys = y*2^25 (a multiple of y, so fmod(x, y) == fmod(x mod ys, y)) to
avoid overflowing when scaling x up. A subnormal x with a normal y is
unaffected: the core takes its ex < ey path and returns x directly.

On flush-to-zero devices __clc_denormals_are_zero_fp32() is true and the
wrapper is skipped, preserving the previous behavior.

Verified against glibc with a scalar model over 40M random pairs with
heavy subnormal coverage (x-sub, y-sub, both): zero mismatches. On
i.MX95 Mali-G310 (Panfrost/Rusticl) advertising CL_FP_DENORM, every
previously failing subnormal point now returns the correct result.

Signed-off-by: jiajia Qian <jiajia.qian at nxp.com>
---
 libclc/clc/lib/generic/math/clc_fmod.cl | 75 ++++++++++++++++++++-----
 1 file changed, 61 insertions(+), 14 deletions(-)

diff --git a/libclc/clc/lib/generic/math/clc_fmod.cl b/libclc/clc/lib/generic/math/clc_fmod.cl
index 7f60b403b53e6..eabfd0ac05204 100644
--- a/libclc/clc/lib/generic/math/clc_fmod.cl
+++ b/libclc/clc/lib/generic/math/clc_fmod.cl
@@ -6,17 +6,23 @@
 //
 //===----------------------------------------------------------------------===//
 
-#include <clc/clc_convert.h>
-#include <clc/integer/clc_clz.h>
-#include <clc/internal/clc.h>
-#include <clc/math/clc_floor.h>
-#include <clc/math/clc_fma.h>
-#include <clc/math/clc_ldexp.h>
-#include <clc/math/clc_trunc.h>
-#include <clc/math/math.h>
-#include <clc/shared/clc_max.h>
-
-_CLC_DEF _CLC_OVERLOAD float __clc_fmod(float x, float y) {
+#include "clc/clc_convert.h"
+#include "clc/integer/clc_clz.h"
+#include "clc/internal/clc.h"
+#include "clc/math/clc_floor.h"
+#include "clc/math/clc_fma.h"
+#include "clc/math/clc_ldexp.h"
+#include "clc/math/clc_subnormal_config.h"
+#include "clc/math/clc_trunc.h"
+#include "clc/math/math.h"
+#include "clc/shared/clc_max.h"
+
+// Core fmod reduction. Assumes x and y are finite, non-zero and normal (the
+// exponent field is reconstructed from ex/ey and the result is rescaled with
+// as_float(ey << 23), which is only correct for normal operands and would
+// flush a subnormal result to zero). Callers must normalize subnormal inputs
+// before calling and handle NaN/Inf/zero edge cases.
+_CLC_DEF _CLC_OVERLOAD float __clc_fmod_impl(float x, float y) {
   int ux = __clc_as_int(x);
   int ax = ux & EXSIGNBIT_SP32;
   float xa = __clc_as_float(ax);
@@ -54,10 +60,51 @@ _CLC_DEF _CLC_OVERLOAD float __clc_fmod(float x, float y) {
   c = ax == ay;
   xr = c ? 0.0f : xr;
 
-  xr = __clc_as_float(sx ^ __clc_as_int(xr));
+  return __clc_as_float(sx ^ __clc_as_int(xr));
+}
+
+_CLC_DEF _CLC_OVERLOAD float __clc_fmod(float x, float y) {
+  int ux = __clc_as_int(x);
+  int ax = ux & EXSIGNBIT_SP32;
+  int uy = __clc_as_int(y);
+  int ay = uy & EXSIGNBIT_SP32;
+
+  float xr;
+
+  // The core reduction only handles normal operands and would flush a
+  // subnormal result to zero (as_float(ey << 23) == 0 when y is subnormal).
+  // When subnormals are supported, normalize a subnormal y by scaling both
+  // operands up by the same power of two: fmod(a*t, b*t) == t*fmod(a,b), so the
+  // core only sees a normal-sized modulus, then scale the result back down.
+  // Scaling by 2^25 turns the largest subnormal normal.
+  bool x_sub = ax != 0 && ax < PINFBITPATT_SP32 &&
+               __clc_as_float(ax) < 0x1.0p-126f;
+  bool y_sub = ay != 0 && ay < PINFBITPATT_SP32 &&
+               __clc_as_float(ay) < 0x1.0p-126f;
+
+  if (!__clc_denormals_are_zero_fp32() && y_sub) {
+    float ys = __clc_ldexp(y, 25);
+    if (x_sub) {
+      // Both subnormal: scaling both up cannot overflow.
+      xr = __clc_ldexp(__clc_fmod_impl(__clc_ldexp(x, 25), ys), -25);
+    } else {
+      // x normal, y subnormal: first reduce x modulo the scaled-up ys = y*2^25
+      // (a multiple of y) with the core; the residue r is in [0, |ys|), i.e.
+      // normal-sized, and satisfies fmod(x, y) == fmod(r, y). Then compute
+      // fmod(r, y) via the scaling identity fmod(r, y) == fmod(r*2^25, ys)/2^25
+      // so the core again only sees normal operands, and scale the result back.
+      float r = __clc_fmod_impl(x, ys);
+      xr = __clc_ldexp(__clc_fmod_impl(__clc_ldexp(r, 25), ys), -25);
+    }
+  } else {
+    // Normal operands, or x subnormal with y normal (the core takes the
+    // ex < ey path and uses the subnormal x directly, without the
+    // as_float(ey << 23) rescale that only applies to the larger operand).
+    xr = __clc_fmod_impl(x, y);
+  }
 
-  c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 |
-      ay == 0;
+  int c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 |
+          ax == PINFBITPATT_SP32 | ay == 0;
   xr = c ? __clc_as_float(QNANBITPATT_SP32) : xr;
 
   return xr;