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 <[email protected]> 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 <[email protected]> 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 _______________________________________________ llvm-branch-commits mailing list [email protected] https://lists.llvm.org/cgi-bin/mailman/listinfo/llvm-branch-commits
