[llvm-branch-commits] [libclc] libclc: Update remainder (PR #187999)
Matt Arsenault via llvm-branch-commits
llvm-branch-commits at lists.llvm.org
Mon Mar 23 02:05:32 PDT 2026
https://github.com/arsenm created https://github.com/llvm/llvm-project/pull/187999
libclc: Update remainder
Previously this was failing conformance without -cl-denorms-are-zero
in the float case, and always failing in the double case.
libclc: Implement remainder with remquo
This fixes conformance failures for double and
without -cl-denorms-are-zero. Optimizations are
able to eliminate the unusued quo handling without
duplicating most of the code.
>From f8c734a25a3549909392af52a89432fe763f3034 Mon Sep 17 00:00:00 2001
From: Matt Arsenault <Matthew.Arsenault at amd.com>
Date: Sun, 22 Mar 2026 10:20:25 +0100
Subject: [PATCH 1/2] libclc: Update remainder
Previously this was failing conformance without -cl-denorms-are-zero
in the float case, and always failing in the double case.
---
libclc/clc/lib/generic/math/clc_remainder.cl | 227 ++----------------
libclc/clc/lib/generic/math/clc_remainder.inc | 171 +++++++++++++
2 files changed, 187 insertions(+), 211 deletions(-)
create mode 100644 libclc/clc/lib/generic/math/clc_remainder.inc
diff --git a/libclc/clc/lib/generic/math/clc_remainder.cl b/libclc/clc/lib/generic/math/clc_remainder.cl
index f1dba87ee5b43..f74d4835d4f2f 100644
--- a/libclc/clc/lib/generic/math/clc_remainder.cl
+++ b/libclc/clc/lib/generic/math/clc_remainder.cl
@@ -6,227 +6,32 @@
//
//===----------------------------------------------------------------------===//
+#include "clc/math/clc_remainder.h"
+
#include "clc/clc_convert.h"
#include "clc/integer/clc_clz.h"
-#include "clc/internal/clc.h"
+
+#include "clc/float/definitions.h"
+#include "clc/math/clc_copysign.h"
+#include "clc/math/clc_fabs.h"
#include "clc/math/clc_floor.h"
+#include "clc/math/clc_flush_if_daz.h"
#include "clc/math/clc_fma.h"
+#include "clc/math/clc_frexp.h"
#include "clc/math/clc_ldexp.h"
-#include "clc/math/clc_remainder.h"
+#include "clc/math/clc_mad.h"
+#include "clc/math/clc_recip_fast.h"
+#include "clc/math/clc_rint.h"
#include "clc/math/clc_trunc.h"
#include "clc/math/math.h"
-#include "clc/shared/clc_max.h"
-
-_CLC_DEF _CLC_OVERLOAD float __clc_remainder(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;
-
- uint q = 0;
-
- while (k > 0) {
- c = xr >= yr;
- q = (q << 1) | c;
- xr -= c ? yr : 0.0f;
- xr += xr;
- --k;
- }
-
- c = xr > yr;
- q = (q << 1) | c;
- xr -= c ? yr : 0.0f;
-
- int lt = ex < ey;
-
- q = lt ? 0 : q;
- xr = lt ? xa : xr;
- yr = lt ? ya : yr;
-
- c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & ((q & 0x1) == 0x1));
- xr -= c ? yr : 0.0f;
- q += c;
-
- float s = __clc_as_float(ey << EXPSHIFTBITS_SP32);
- xr *= lt ? 1.0f : s;
+#include "clc/relational/clc_isfinite.h"
+#include "clc/relational/clc_isnan.h"
+#include "clc/relational/clc_signbit.h"
- 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;
-}
-
-#define __CLC_FLOAT_ONLY
-#define __CLC_FUNCTION __clc_remainder
-#define __CLC_BODY "clc/shared/binary_def_scalarize.inc"
+#define __CLC_BODY "clc_remainder.inc"
#include "clc/math/gentype.inc"
-#undef __CLC_FUNCTION
-
-#ifdef cl_khr_fp64
-
-#pragma OPENCL EXTENSION cl_khr_fp64 : enable
-
-_CLC_DEF _CLC_OVERLOAD double __clc_remainder(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;
-
- int qsgn = ((ux ^ uy) & SIGNBIT_DP64) == 0UL ? 1 : -1;
-
- // 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;
- }
-
- // One more time
- // Variable todd says whether the integer t is odd or not
- t = __clc_floor(dx / w);
- long lt = (long)t;
- int todd = lt & 1;
-
- p = w * t;
- pp = __clc_fma(w, t, -p);
- v = dx - p;
- dx = v + (((dx - v) - p) - pp);
- i = dx < 0.0;
- todd ^= i;
- dx += i ? w : 0.0;
-
- // At this point, dx lies in the range [0,dy)
-
- // For the fmod function, we're done apart from setting the correct sign.
- //
- // For the remainder function, we need to adjust dx
- // so that it lies in the range (-y/2, y/2] by carefully
- // subtracting w (== dy == y) if necessary. The rigmarole
- // with todd is to get the correct sign of the result
- // when x/y lies exactly half way between two integers,
- // when we need to choose the even integer.
-
- int al = (2.0 * dx > w) | (todd & (2.0 * dx == w));
- double dxl = dx - (al ? w : 0.0);
-
- int ag = (dx > 0.5 * w) | (todd & (dx == 0.5 * w));
- double dxg = dx - (ag ? w : 0.0);
-
- dx = dy < 0x1.0p+1022 ? dxl : dxg;
-
- 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;
-
- c &= (yexp<1023 & 2.0 * dx> dy) | (dx > 0.5 * dy);
- // we could use a conversion here instead since qsgn = +-1
- p = qsgn == 1 ? -1.0 : 1.0;
- t = __clc_fma(y, p, x);
- ret = c ? t : 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;
-
- return ret;
-}
-
-#define __CLC_DOUBLE_ONLY
#define __CLC_FUNCTION __clc_remainder
-#define __CLC_BODY "clc/shared/binary_def_scalarize.inc"
+#define __CLC_BODY "clc/shared/binary_def_scalarize_loop.inc"
#include "clc/math/gentype.inc"
#undef __CLC_FUNCTION
-
-#endif
-
-#ifdef cl_khr_fp16
-
-#pragma OPENCL EXTENSION cl_khr_fp16 : enable
-
-// Forward the half version of this builtin onto the float one
-#define __CLC_HALF_ONLY
-#define __CLC_FUNCTION __clc_remainder
-#define __CLC_BODY "clc/math/binary_def_via_fp32.inc"
-#include "clc/math/gentype.inc"
-
-#endif
diff --git a/libclc/clc/lib/generic/math/clc_remainder.inc b/libclc/clc/lib/generic/math/clc_remainder.inc
new file mode 100644
index 0000000000000..eba99750e082e
--- /dev/null
+++ b/libclc/clc/lib/generic/math/clc_remainder.inc
@@ -0,0 +1,171 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifdef __CLC_SCALAR
+
+#if __CLC_FPSIZE == 32
+#define __CLC_REMAINDER_EVAL_TYPE __CLC_FLOATN
+#define __CLC_CONVERT_REMAINDER_EVAL_TYPE __CLC_CONVERT_FLOATN
+#define __CLC_S_EVAL_TYPE __CLC_INTN
+#define __CLC_CONVERT_S_EVAL_TYPE __CLC_CONVERT_INTN
+#elif __CLC_FPSIZE == 64
+#define __CLC_REMAINDER_EVAL_TYPE __CLC_DOUBLEN
+#define __CLC_CONVERT_REMAINDER_EVAL_TYPE __CLC_CONVERT_DOUBLEN
+#define __CLC_S_EVAL_TYPE __CLC_LONGN
+#define __CLC_CONVERT_S_EVAL_TYPE __CLC_CONVERT_LONGN
+#elif __CLC_FPSIZE == 16
+#define __CLC_REMAINDER_EVAL_TYPE __CLC_FLOATN
+#define __CLC_CONVERT_REMAINDER_EVAL_TYPE __CLC_CONVERT_FLOATN
+#define __CLC_S_EVAL_TYPE __CLC_INTN
+#define __CLC_CONVERT_S_EVAL_TYPE __CLC_CONVERT_INTN
+#endif
+
+// The arguments must only be variable names
+#define FULL_MUL(A, B, CHI, CLO) \
+ do { \
+ __CLC_GENTYPE __ha = \
+ __CLC_AS_GENTYPE(__CLC_AS_U_GENTYPE(A) & 0xfffff000U); \
+ __CLC_GENTYPE __ta = A - __ha; \
+ __CLC_GENTYPE __hb = \
+ __CLC_AS_GENTYPE(__CLC_AS_U_GENTYPE(B) & 0xfffff000U); \
+ __CLC_GENTYPE __tb = B - __hb; \
+ CHI = A * B; \
+ CLO = __clc_mad( \
+ __ta, __tb, \
+ __clc_mad(__ta, __hb, \
+ __clc_mad(__ha, __tb, __clc_mad(__ha, __hb, -CHI)))); \
+ } while (0)
+
+#if __CLC_FPSIZE != 16
+static _CLC_OVERLOAD _CLC_CONST __CLC_REMAINDER_EVAL_TYPE
+fnma(__CLC_REMAINDER_EVAL_TYPE a, __CLC_REMAINDER_EVAL_TYPE b,
+ __CLC_REMAINDER_EVAL_TYPE c) {
+#if __CLC_FPSIZE == 32
+ if (false && !HAVE_FAST_FMA32()) {
+ __CLC_GENTYPE h, t;
+ FULL_MUL(a, b, h, t);
+ __CLC_GENTYPE d = c - h;
+ d = (((c - d) - h) - t) + d;
+ return d;
+ }
+#endif
+
+ return __clc_fma(-a, b, c);
+}
+#endif
+
+_CLC_DEF _CLC_OVERLOAD _CLC_CONST __CLC_GENTYPE
+__clc_remainder(__CLC_GENTYPE x, __CLC_GENTYPE y) {
+ // How many bits of the quotient per iteration
+
+#if __CLC_FPSIZE == 32
+ const __CLC_INTN bits = 12;
+ const __CLC_GENTYPE max_exp = 0x1.0p+127f;
+#elif __CLC_FPSIZE == 64
+ const __CLC_INTN bits = 26;
+ const __CLC_GENTYPE max_exp = 0x1.0p+1023;
+#elif __CLC_FPSIZE == 16
+ const __CLC_INTN bits = 11;
+ const __CLC_GENTYPE max_exp = 0x1.0p+15h;
+#endif
+
+ __CLC_REMAINDER_EVAL_TYPE ax =
+ __CLC_CONVERT_REMAINDER_EVAL_TYPE(__clc_fabs(x));
+ __CLC_REMAINDER_EVAL_TYPE ay =
+ __CLC_CONVERT_REMAINDER_EVAL_TYPE(__clc_fabs(y));
+
+ __CLC_GENTYPE ret;
+
+ if (ax > ay) {
+ __CLC_INTN ex, ey;
+
+ __CLC_REMAINDER_EVAL_TYPE mx = __clc_frexp(ax, &ex);
+ --ex;
+
+ __CLC_REMAINDER_EVAL_TYPE my = __clc_frexp(ay, &ey);
+ --ey;
+
+ ax = __clc_ldexp(mx, bits);
+ ay = __clc_ldexp(my, 1);
+
+ __CLC_INTN nb = ex - ey;
+ __CLC_REMAINDER_EVAL_TYPE ayinv = __clc_recip_fast(ay);
+
+ __CLC_INTN qacc = 0;
+
+ while (nb > bits) {
+ __CLC_REMAINDER_EVAL_TYPE q = __clc_rint(ax * ayinv);
+
+#if __CLC_FPSIZE == 16
+ ax = __clc_mad(-q, ay, ax);
+#else
+ ax = fnma(q, ay, ax);
+#endif
+ __CLC_S_GENTYPE clt = ax < (__CLC_REMAINDER_EVAL_TYPE)0.0;
+ __CLC_REMAINDER_EVAL_TYPE axp = ax + ay;
+ ax = clt ? axp : ax;
+ ax = __clc_ldexp(ax, bits);
+ nb -= bits;
+ }
+
+ ax = __clc_ldexp(ax, nb - bits + 1);
+
+ // Final iteration
+ {
+ __CLC_REMAINDER_EVAL_TYPE q = __clc_rint(ax * ayinv);
+#if __CLC_FPSIZE == 16
+ ax = __clc_mad(-q, ay, ax);
+#else
+ ax = fnma(q, ay, ax);
+#endif
+ __CLC_S_GENTYPE clt = ax < (__CLC_REMAINDER_EVAL_TYPE)0.0;
+ __CLC_REMAINDER_EVAL_TYPE axp = ax + ay;
+ ax = clt ? axp : ax;
+ __CLC_INTN iq = __CLC_CONVERT_INTN(q);
+ iq -= __CLC_CONVERT_INTN(clt) ? 1 : 0;
+ qacc = iq;
+ }
+
+ // Adjust ax so that it is the range (-y/2, y/2]
+ // We need to choose the even integer when x/y is midway between two
+ // integers
+ __CLC_S_EVAL_TYPE aq = ((__CLC_REMAINDER_EVAL_TYPE)2.0 * ax > ay) |
+ (__CLC_CONVERT_S_EVAL_TYPE(qacc & 0x1) &
+ ((__CLC_REMAINDER_EVAL_TYPE)2.0 * ax == ay));
+ ax = ax - (aq ? ay : (__CLC_REMAINDER_EVAL_TYPE)0.0);
+
+ ax = __clc_ldexp(ax, ey);
+
+ ret = __clc_signbit(x) ? -ax : ax;
+ } else {
+ __CLC_S_EVAL_TYPE c = (ax > (__CLC_REMAINDER_EVAL_TYPE)0.5 * ay);
+ if (__CLC_FPSIZE != 16)
+ c |= (ay < max_exp && (__CLC_REMAINDER_EVAL_TYPE)2.0 * ax > ay);
+
+ __CLC_CHARN qsgn = __CLC_CONVERT_CHARN(__clc_signbit(x) == __clc_signbit(y))
+ ? (__CLC_CHARN)1
+ : (__CLC_CHARN)-1;
+
+ __CLC_GENTYPE t = __clc_mad(y, -__CLC_CONVERT_GENTYPE(qsgn), x);
+ ret = c ? t : __clc_flush_if_daz(x);
+
+ __CLC_GENTYPE zero = __clc_copysign(__CLC_FP_LIT(0.0), x);
+ ret = ax == ay ? zero : ret;
+ }
+
+ ret = y == __CLC_FP_LIT(0.0) ? __CLC_GENTYPE_NAN : ret;
+
+ return !__clc_isnan(y) && __clc_isfinite(x) ? ret : __CLC_GENTYPE_NAN;
+}
+
+#undef __CLC_REMAINDER_EVAL_TYPE
+#undef __CLC_CONVERT_REMAINDER_EVAL_TYPE
+#undef __CLC_S_EVAL_TYPE
+#undef __CLC_CONVERT_S_EVAL_TYPE
+
+#endif // __CLC_SCALAR
>From b6deaa6eace2bed85a502eee76caef24fc57e6ce Mon Sep 17 00:00:00 2001
From: Matt Arsenault <Matthew.Arsenault at amd.com>
Date: Sun, 22 Mar 2026 19:21:08 +0100
Subject: [PATCH 2/2] libclc: Implement remainder with remquo
This fixes conformance failures for double and
without -cl-denorms-are-zero. Optimizations are
able to eliminate the unusued quo handling without
duplicating most of the code.
---
libclc/clc/lib/generic/math/clc_remainder.cl | 26 +--
libclc/clc/lib/generic/math/clc_remainder.inc | 162 +-----------------
2 files changed, 3 insertions(+), 185 deletions(-)
diff --git a/libclc/clc/lib/generic/math/clc_remainder.cl b/libclc/clc/lib/generic/math/clc_remainder.cl
index f74d4835d4f2f..d3979fbac3ffd 100644
--- a/libclc/clc/lib/generic/math/clc_remainder.cl
+++ b/libclc/clc/lib/generic/math/clc_remainder.cl
@@ -7,31 +7,7 @@
//===----------------------------------------------------------------------===//
#include "clc/math/clc_remainder.h"
-
-#include "clc/clc_convert.h"
-#include "clc/integer/clc_clz.h"
-
-#include "clc/float/definitions.h"
-#include "clc/math/clc_copysign.h"
-#include "clc/math/clc_fabs.h"
-#include "clc/math/clc_floor.h"
-#include "clc/math/clc_flush_if_daz.h"
-#include "clc/math/clc_fma.h"
-#include "clc/math/clc_frexp.h"
-#include "clc/math/clc_ldexp.h"
-#include "clc/math/clc_mad.h"
-#include "clc/math/clc_recip_fast.h"
-#include "clc/math/clc_rint.h"
-#include "clc/math/clc_trunc.h"
-#include "clc/math/math.h"
-#include "clc/relational/clc_isfinite.h"
-#include "clc/relational/clc_isnan.h"
-#include "clc/relational/clc_signbit.h"
+#include "clc/math/clc_remquo.h"
#define __CLC_BODY "clc_remainder.inc"
#include "clc/math/gentype.inc"
-
-#define __CLC_FUNCTION __clc_remainder
-#define __CLC_BODY "clc/shared/binary_def_scalarize_loop.inc"
-#include "clc/math/gentype.inc"
-#undef __CLC_FUNCTION
diff --git a/libclc/clc/lib/generic/math/clc_remainder.inc b/libclc/clc/lib/generic/math/clc_remainder.inc
index eba99750e082e..9f2e0ba1539a6 100644
--- a/libclc/clc/lib/generic/math/clc_remainder.inc
+++ b/libclc/clc/lib/generic/math/clc_remainder.inc
@@ -6,166 +6,8 @@
//
//===----------------------------------------------------------------------===//
-#ifdef __CLC_SCALAR
-
-#if __CLC_FPSIZE == 32
-#define __CLC_REMAINDER_EVAL_TYPE __CLC_FLOATN
-#define __CLC_CONVERT_REMAINDER_EVAL_TYPE __CLC_CONVERT_FLOATN
-#define __CLC_S_EVAL_TYPE __CLC_INTN
-#define __CLC_CONVERT_S_EVAL_TYPE __CLC_CONVERT_INTN
-#elif __CLC_FPSIZE == 64
-#define __CLC_REMAINDER_EVAL_TYPE __CLC_DOUBLEN
-#define __CLC_CONVERT_REMAINDER_EVAL_TYPE __CLC_CONVERT_DOUBLEN
-#define __CLC_S_EVAL_TYPE __CLC_LONGN
-#define __CLC_CONVERT_S_EVAL_TYPE __CLC_CONVERT_LONGN
-#elif __CLC_FPSIZE == 16
-#define __CLC_REMAINDER_EVAL_TYPE __CLC_FLOATN
-#define __CLC_CONVERT_REMAINDER_EVAL_TYPE __CLC_CONVERT_FLOATN
-#define __CLC_S_EVAL_TYPE __CLC_INTN
-#define __CLC_CONVERT_S_EVAL_TYPE __CLC_CONVERT_INTN
-#endif
-
-// The arguments must only be variable names
-#define FULL_MUL(A, B, CHI, CLO) \
- do { \
- __CLC_GENTYPE __ha = \
- __CLC_AS_GENTYPE(__CLC_AS_U_GENTYPE(A) & 0xfffff000U); \
- __CLC_GENTYPE __ta = A - __ha; \
- __CLC_GENTYPE __hb = \
- __CLC_AS_GENTYPE(__CLC_AS_U_GENTYPE(B) & 0xfffff000U); \
- __CLC_GENTYPE __tb = B - __hb; \
- CHI = A * B; \
- CLO = __clc_mad( \
- __ta, __tb, \
- __clc_mad(__ta, __hb, \
- __clc_mad(__ha, __tb, __clc_mad(__ha, __hb, -CHI)))); \
- } while (0)
-
-#if __CLC_FPSIZE != 16
-static _CLC_OVERLOAD _CLC_CONST __CLC_REMAINDER_EVAL_TYPE
-fnma(__CLC_REMAINDER_EVAL_TYPE a, __CLC_REMAINDER_EVAL_TYPE b,
- __CLC_REMAINDER_EVAL_TYPE c) {
-#if __CLC_FPSIZE == 32
- if (false && !HAVE_FAST_FMA32()) {
- __CLC_GENTYPE h, t;
- FULL_MUL(a, b, h, t);
- __CLC_GENTYPE d = c - h;
- d = (((c - d) - h) - t) + d;
- return d;
- }
-#endif
-
- return __clc_fma(-a, b, c);
-}
-#endif
-
_CLC_DEF _CLC_OVERLOAD _CLC_CONST __CLC_GENTYPE
__clc_remainder(__CLC_GENTYPE x, __CLC_GENTYPE y) {
- // How many bits of the quotient per iteration
-
-#if __CLC_FPSIZE == 32
- const __CLC_INTN bits = 12;
- const __CLC_GENTYPE max_exp = 0x1.0p+127f;
-#elif __CLC_FPSIZE == 64
- const __CLC_INTN bits = 26;
- const __CLC_GENTYPE max_exp = 0x1.0p+1023;
-#elif __CLC_FPSIZE == 16
- const __CLC_INTN bits = 11;
- const __CLC_GENTYPE max_exp = 0x1.0p+15h;
-#endif
-
- __CLC_REMAINDER_EVAL_TYPE ax =
- __CLC_CONVERT_REMAINDER_EVAL_TYPE(__clc_fabs(x));
- __CLC_REMAINDER_EVAL_TYPE ay =
- __CLC_CONVERT_REMAINDER_EVAL_TYPE(__clc_fabs(y));
-
- __CLC_GENTYPE ret;
-
- if (ax > ay) {
- __CLC_INTN ex, ey;
-
- __CLC_REMAINDER_EVAL_TYPE mx = __clc_frexp(ax, &ex);
- --ex;
-
- __CLC_REMAINDER_EVAL_TYPE my = __clc_frexp(ay, &ey);
- --ey;
-
- ax = __clc_ldexp(mx, bits);
- ay = __clc_ldexp(my, 1);
-
- __CLC_INTN nb = ex - ey;
- __CLC_REMAINDER_EVAL_TYPE ayinv = __clc_recip_fast(ay);
-
- __CLC_INTN qacc = 0;
-
- while (nb > bits) {
- __CLC_REMAINDER_EVAL_TYPE q = __clc_rint(ax * ayinv);
-
-#if __CLC_FPSIZE == 16
- ax = __clc_mad(-q, ay, ax);
-#else
- ax = fnma(q, ay, ax);
-#endif
- __CLC_S_GENTYPE clt = ax < (__CLC_REMAINDER_EVAL_TYPE)0.0;
- __CLC_REMAINDER_EVAL_TYPE axp = ax + ay;
- ax = clt ? axp : ax;
- ax = __clc_ldexp(ax, bits);
- nb -= bits;
- }
-
- ax = __clc_ldexp(ax, nb - bits + 1);
-
- // Final iteration
- {
- __CLC_REMAINDER_EVAL_TYPE q = __clc_rint(ax * ayinv);
-#if __CLC_FPSIZE == 16
- ax = __clc_mad(-q, ay, ax);
-#else
- ax = fnma(q, ay, ax);
-#endif
- __CLC_S_GENTYPE clt = ax < (__CLC_REMAINDER_EVAL_TYPE)0.0;
- __CLC_REMAINDER_EVAL_TYPE axp = ax + ay;
- ax = clt ? axp : ax;
- __CLC_INTN iq = __CLC_CONVERT_INTN(q);
- iq -= __CLC_CONVERT_INTN(clt) ? 1 : 0;
- qacc = iq;
- }
-
- // Adjust ax so that it is the range (-y/2, y/2]
- // We need to choose the even integer when x/y is midway between two
- // integers
- __CLC_S_EVAL_TYPE aq = ((__CLC_REMAINDER_EVAL_TYPE)2.0 * ax > ay) |
- (__CLC_CONVERT_S_EVAL_TYPE(qacc & 0x1) &
- ((__CLC_REMAINDER_EVAL_TYPE)2.0 * ax == ay));
- ax = ax - (aq ? ay : (__CLC_REMAINDER_EVAL_TYPE)0.0);
-
- ax = __clc_ldexp(ax, ey);
-
- ret = __clc_signbit(x) ? -ax : ax;
- } else {
- __CLC_S_EVAL_TYPE c = (ax > (__CLC_REMAINDER_EVAL_TYPE)0.5 * ay);
- if (__CLC_FPSIZE != 16)
- c |= (ay < max_exp && (__CLC_REMAINDER_EVAL_TYPE)2.0 * ax > ay);
-
- __CLC_CHARN qsgn = __CLC_CONVERT_CHARN(__clc_signbit(x) == __clc_signbit(y))
- ? (__CLC_CHARN)1
- : (__CLC_CHARN)-1;
-
- __CLC_GENTYPE t = __clc_mad(y, -__CLC_CONVERT_GENTYPE(qsgn), x);
- ret = c ? t : __clc_flush_if_daz(x);
-
- __CLC_GENTYPE zero = __clc_copysign(__CLC_FP_LIT(0.0), x);
- ret = ax == ay ? zero : ret;
- }
-
- ret = y == __CLC_FP_LIT(0.0) ? __CLC_GENTYPE_NAN : ret;
-
- return !__clc_isnan(y) && __clc_isfinite(x) ? ret : __CLC_GENTYPE_NAN;
+ __CLC_INTN unused_quo;
+ return __clc_remquo(x, y, &unused_quo);
}
-
-#undef __CLC_REMAINDER_EVAL_TYPE
-#undef __CLC_CONVERT_REMAINDER_EVAL_TYPE
-#undef __CLC_S_EVAL_TYPE
-#undef __CLC_CONVERT_S_EVAL_TYPE
-
-#endif // __CLC_SCALAR
More information about the llvm-branch-commits
mailing list