[libclc] ab7245d - libclc: Update fmod implementations (#222369)
via cfe-commits
cfe-commits at lists.llvm.org
Thu Sep 10 00:56:21 PDT 2026
Author: Matt Arsenault
Date: 2026-09-10T09:56:16+02:00
New Revision: ab7245d75cc97db490d73037b5f021a561381d0d
URL: https://github.com/llvm/llvm-project/commit/ab7245d75cc97db490d73037b5f021a561381d0d
DIFF: https://github.com/llvm/llvm-project/commit/ab7245d75cc97db490d73037b5f021a561381d0d.diff
LOG: 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 at anthropic.com>
Added:
Modified:
libclc/clc/lib/generic/math/clc_fmod.cl
Removed:
################################################################################
diff --git a/libclc/clc/lib/generic/math/clc_fmod.cl b/libclc/clc/lib/generic/math/clc_fmod.cl
index 7f60b403b53e6..699d95485b71e 100644
--- a/libclc/clc/lib/generic/math/clc_fmod.cl
+++ b/libclc/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;
-
- int uy = __clc_as_int(y);
- int ay = uy & EXSIGNBIT_SP32;
- float ya = __clc_as_float(ay);
- int ey = ay >> EXPSHIFTBITS_SP32;
-
- float xr = __clc_as_float(0x3f800000 | (ax & 0x007fffff));
- float yr = __clc_as_float(0x3f800000 | (ay & 0x007fffff));
- int c;
- int k = ex - ey;
-
- while (k > 0) {
- c = xr >= yr;
- xr -= c ? yr : 0.0f;
- xr += xr;
- --k;
+ // How many bits of the quotient to resolve per iteration.
+ const int bits = 12;
+
+ float ax = __clc_fabs(x);
+ float ay = __clc_fabs(y);
+
+ float ret;
+
+ 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 @@ _CLC_DEF _CLC_OVERLOAD float __clc_fmod(float x, float y) {
#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;
-
- 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;
-
- // First assume |x| > |y|
-
- // 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;
-
- // 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;
-
- for (i = 0; i < ntimes; i++) {
- // Compute integral multiplier
- t = __clc_trunc(dx / w);
-
- // Compute w * t in quad precision
- p = w * t;
- pp = __clc_fma(w, t, -p);
-
- // Subtract w * t from dx
- v = dx - p;
- dx = v + (((dx - v) - p) - pp);
-
- // If t was one too large, dx will be negative. Add back one w.
- dx += dx < 0.0 ? w : 0.0;
-
- // Scale w down by 2^(-53) for the next iteration
- w *= scale;
+ // How many bits of the quotient to resolve per iteration.
+ const int bits = 26;
+
+ double ax = __clc_fabs(x);
+ double ay = __clc_fabs(y);
+
+ double ret;
+
+ if (ax > ay) {
+ int ex, ey;
+ double mx = __clc_frexp(ax, &ex);
+ --ex;
+ double my = __clc_frexp(ay, &ey);
+ --ey;
+
+ ax = __clc_ldexp(mx, bits);
+ ay = __clc_ldexp(my, 1);
+
+ int nb = ex - ey;
+ double ayinv = 1.0 / ay;
+
+ 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;
+ }
+
+ ax = __clc_ldexp(ax, nb - bits + 1);
+
+ // 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;
+
+ 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;
}
More information about the cfe-commits
mailing list