Google GitGoogle Git
Sign in
llvm/llvm-project/libclc/refs/heads/master^!/.
commit
5932501d8a1543d8c4796a04e65596fc5f81f621
[log]
author
Matt Arsenault <Matthew.Arsenault@amd.com>
Thu Sep 10 09:56:16 2026 +0200
committer
Copybara-Service <copybara-worker@google.com>
Thu Sep 10 01:03:18 2026 -0700
tree
478e3651a46109bff1b1df160d5c375935d9501f
parent
12e4086b6dc0aa77e121467b58a5a02b2f3a532d [diff]
libclc: Update fmod implementations (#222369)

This was originally ported from rocm device libs in
93af966747b59d37c57312a0c0242151076c072b. Merge in more
recent changes. This should also approximately match the default
expansion in ExpandIRInsts

Co-authored-by: Claude <noreply@anthropic.com>
GitOrigin-RevId: ab7245d75cc97db490d73037b5f021a561381d0d

diff --git a/clc/lib/generic/math/clc_fmod.cl b/clc/lib/generic/math/clc_fmod.cl
index 7f60b40..699d954 100644
--- a/clc/lib/generic/math/clc_fmod.cl
+++ b/clc/lib/generic/math/clc_fmod.cl

@@ -7,60 +7,73 @@
 //===----------------------------------------------------------------------===//
 
 #include <clc/clc_convert.h>
-#include <clc/integer/clc_clz.h>
+#include <clc/float/definitions.h>
 #include <clc/internal/clc.h>
-#include <clc/math/clc_floor.h>
+#include <clc/math/clc_copysign.h>
+#include <clc/math/clc_fabs.h>
 #include <clc/math/clc_fma.h>
+#include <clc/math/clc_frexp.h>
 #include <clc/math/clc_ldexp.h>
-#include <clc/math/clc_trunc.h>
+#include <clc/math/clc_recip_fast.h>
+#include <clc/math/clc_rint.h>
 #include <clc/math/math.h>
-#include <clc/shared/clc_max.h>
+#include <clc/relational/clc_isfinite.h>
+#include <clc/relational/clc_isnan.h>
 
 _CLC_DEF _CLC_OVERLOAD float __clc_fmod(float x, float y) {
-  int ux = __clc_as_int(x);
-  int ax = ux & EXSIGNBIT_SP32;
-  float xa = __clc_as_float(ax);
-  int sx = ux ^ ax;
-  int ex = ax >> EXPSHIFTBITS_SP32;
+  // How many bits of the quotient to resolve per iteration.
+  const int bits = 12;
 
-  int uy = __clc_as_int(y);
-  int ay = uy & EXSIGNBIT_SP32;
-  float ya = __clc_as_float(ay);
-  int ey = ay >> EXPSHIFTBITS_SP32;
+  float ax = __clc_fabs(x);
+  float ay = __clc_fabs(y);
 
-  float xr = __clc_as_float(0x3f800000 | (ax & 0x007fffff));
-  float yr = __clc_as_float(0x3f800000 | (ay & 0x007fffff));
-  int c;
-  int k = ex - ey;
+  float ret;
 
-  while (k > 0) {
-    c = xr >= yr;
-    xr -= c ? yr : 0.0f;
-    xr += xr;
-    --k;
+  if (ax > ay) {
+    int ex, ey;
+    float mx = __clc_frexp(ax, &ex);
+    --ex;
+    float my = __clc_frexp(ay, &ey);
+    --ey;
+
+    ax = __clc_ldexp(mx, bits);
+    ay = __clc_ldexp(my, 1);
+
+    int nb = ex - ey;
+    float ayinv = __clc_recip_fast(ay);
+
+    while (nb > bits) {
+      float q = __clc_rint(ax * ayinv);
+      ax = __clc_fma(-q, ay, ax);
+      int clt = ax < 0.0f;
+      float axp = ax + ay;
+      ax = clt ? axp : ax;
+      ax = __clc_ldexp(ax, bits);
+      nb -= bits;
+    }
+
+    ax = __clc_ldexp(ax, nb - bits + 1);
+
+    // Final iteration.
+    float q = __clc_rint(ax * ayinv);
+    ax = __clc_fma(-q, ay, ax);
+    int clt = ax < 0.0f;
+    float axp = ax + ay;
+    ax = clt ? axp : ax;
+
+    ax = __clc_ldexp(ax, ey);
+    ret = __clc_as_float((__clc_as_int(x) & SIGNBIT_SP32) ^ __clc_as_int(ax));
+  } else {
+    // |x| < |y| returns x; |x| == |y| returns a zero with the sign of x.
+    ret = ax == ay ? __clc_copysign(0.0f, x) : x;
   }
 
-  c = xr >= yr;
-  xr -= c ? yr : 0.0f;
+  // fmod(x, 0) is NaN; fmod(Inf, y) is NaN; fmod(x, NaN)/fmod(NaN, y) is NaN.
+  ret = y == 0.0f ? FLT_NAN : ret;
+  int c = !__clc_isnan(y) && __clc_isfinite(x);
+  ret = c ? ret : FLT_NAN;
 
-  int lt = ex < ey;
-
-  xr = lt ? xa : xr;
-  yr = lt ? ya : yr;
-
-  float s = __clc_as_float(ey << EXPSHIFTBITS_SP32);
-  xr *= lt ? 1.0f : s;
-
-  c = ax == ay;
-  xr = c ? 0.0f : xr;
-
-  xr = __clc_as_float(sx ^ __clc_as_int(xr));
-
-  c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 |
-      ay == 0;
-  xr = c ? __clc_as_float(QNANBITPATT_SP32) : xr;
-
-  return xr;
+  return ret;
 }
 
 #define __CLC_FLOAT_ONLY
@@ -74,99 +87,58 @@
 #pragma OPENCL EXTENSION cl_khr_fp64 : enable
 
 _CLC_DEF _CLC_OVERLOAD double __clc_fmod(double x, double y) {
-  ulong ux = __clc_as_ulong(x);
-  ulong ax = ux & ~SIGNBIT_DP64;
-  ulong xsgn = ux ^ ax;
-  double dx = __clc_as_double(ax);
-  int xexp = __clc_convert_int(ax >> EXPSHIFTBITS_DP64);
-  int xexp1 = 11 - (int)__clc_clz(ax & MANTBITS_DP64);
-  xexp1 = xexp < 1 ? xexp1 : xexp;
+  // How many bits of the quotient to resolve per iteration.
+  const int bits = 26;
 
-  ulong uy = __clc_as_ulong(y);
-  ulong ay = uy & ~SIGNBIT_DP64;
-  double dy = __clc_as_double(ay);
-  int yexp = __clc_convert_int(ay >> EXPSHIFTBITS_DP64);
-  int yexp1 = 11 - (int)__clc_clz(ay & MANTBITS_DP64);
-  yexp1 = yexp < 1 ? yexp1 : yexp;
+  double ax = __clc_fabs(x);
+  double ay = __clc_fabs(y);
 
-  // First assume |x| > |y|
+  double ret;
 
-  // Set ntimes to the number of times we need to do a
-  // partial remainder. If the exponent of x is an exact multiple
-  // of 53 larger than the exponent of y, and the mantissa of x is
-  // less than the mantissa of y, ntimes will be one too large
-  // but it doesn't matter - it just means that we'll go round
-  // the loop below one extra time.
-  int ntimes = __clc_max(0, (xexp1 - yexp1) / 53);
-  double w = __clc_ldexp(dy, ntimes * 53);
-  w = ntimes == 0 ? dy : w;
-  double scale = ntimes == 0 ? 1.0 : 0x1.0p-53;
+  if (ax > ay) {
+    int ex, ey;
+    double mx = __clc_frexp(ax, &ex);
+    --ex;
+    double my = __clc_frexp(ay, &ey);
+    --ey;
 
-  // Each time round the loop we compute a partial remainder.
-  // This is done by subtracting a large multiple of w
-  // from x each time, where w is a scaled up version of y.
-  // The subtraction must be performed exactly in quad
-  // precision, though the result at each stage can
-  // fit exactly in a double precision number.
-  int i;
-  double t, v, p, pp;
+    ax = __clc_ldexp(mx, bits);
+    ay = __clc_ldexp(my, 1);
 
-  for (i = 0; i < ntimes; i++) {
-    // Compute integral multiplier
-    t = __clc_trunc(dx / w);
+    int nb = ex - ey;
+    double ayinv = 1.0 / ay;
 
-    // Compute w * t in quad precision
-    p = w * t;
-    pp = __clc_fma(w, t, -p);
+    while (nb > bits) {
+      double q = __clc_rint(ax * ayinv);
+      ax = __clc_fma(-q, ay, ax);
+      int clt = ax < 0.0;
+      double axp = ax + ay;
+      ax = clt ? axp : ax;
+      ax = __clc_ldexp(ax, bits);
+      nb -= bits;
+    }
 
-    // Subtract w * t from dx
-    v = dx - p;
-    dx = v + (((dx - v) - p) - pp);
+    ax = __clc_ldexp(ax, nb - bits + 1);
 
-    // If t was one too large, dx will be negative. Add back one w.
-    dx += dx < 0.0 ? w : 0.0;
+    // Final iteration.
+    double q = __clc_rint(ax * ayinv);
+    ax = __clc_fma(-q, ay, ax);
+    int clt = ax < 0.0;
+    double axp = ax + ay;
+    ax = clt ? axp : ax;
 
-    // Scale w down by 2^(-53) for the next iteration
-    w *= scale;
+    ax = __clc_ldexp(ax, ey);
+    ret = __clc_as_double((__clc_as_ulong(x) & SIGNBIT_DP64) ^
+                          __clc_as_ulong(ax));
+  } else {
+    // |x| < |y| returns x; |x| == |y| returns a zero with the sign of x.
+    ret = ax == ay ? __clc_copysign(0.0, x) : x;
   }
 
-  // One more time
-  t = __clc_floor(dx / w);
-
-  p = w * t;
-  pp = __clc_fma(w, t, -p);
-  v = dx - p;
-  dx = v + (((dx - v) - p) - pp);
-  i = dx < 0.0;
-  dx += i ? w : 0.0;
-
-  // At this point, dx lies in the range [0,dy)
-  double ret = __clc_as_double(xsgn ^ __clc_as_ulong(dx));
-  dx = __clc_as_double(ax);
-
-  // Now handle |x| == |y|
-  int c = dx == dy;
-  t = __clc_as_double(xsgn);
-  ret = c ? t : ret;
-
-  // Next, handle |x| < |y|
-  c = dx < dy;
-  ret = c ? x : ret;
-
-  // We don't need anything special for |x| == 0
-
-  // |y| is 0
-  c = dy == 0.0;
-  ret = c ? __clc_as_double(QNANBITPATT_DP64) : ret;
-
-  // y is +-Inf, NaN
-  c = yexp > BIASEDEMAX_DP64;
-  t = y == y ? x : y;
-  ret = c ? t : ret;
-
-  // x is +=Inf, NaN
-  c = xexp > BIASEDEMAX_DP64;
-  ret = c ? __clc_as_double(QNANBITPATT_DP64) : ret;
+  // fmod(x, 0) is NaN; fmod(Inf, y) is NaN; fmod(x, NaN)/fmod(NaN, y) is NaN.
+  ret = y == 0.0 ? DBL_NAN : ret;
+  int c = !__clc_isnan(y) && __clc_isfinite(x);
+  ret = c ? ret : DBL_NAN;
 
   return ret;
 }