233 lines · plain
1//===----------------------------------------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8 9#include <clc/clc_convert.h>10#include <clc/integer/clc_clz.h>11#include <clc/internal/clc.h>12#include <clc/math/clc_floor.h>13#include <clc/math/clc_fma.h>14#include <clc/math/clc_ldexp.h>15#include <clc/math/clc_remainder.h>16#include <clc/math/clc_trunc.h>17#include <clc/math/math.h>18#include <clc/shared/clc_max.h>19 20_CLC_DEF _CLC_OVERLOAD float __clc_remainder(float x, float y) {21 int ux = __clc_as_int(x);22 int ax = ux & EXSIGNBIT_SP32;23 float xa = __clc_as_float(ax);24 int sx = ux ^ ax;25 int ex = ax >> EXPSHIFTBITS_SP32;26 27 int uy = __clc_as_int(y);28 int ay = uy & EXSIGNBIT_SP32;29 float ya = __clc_as_float(ay);30 int ey = ay >> EXPSHIFTBITS_SP32;31 32 float xr = __clc_as_float(0x3f800000 | (ax & 0x007fffff));33 float yr = __clc_as_float(0x3f800000 | (ay & 0x007fffff));34 int c;35 int k = ex - ey;36 37 uint q = 0;38 39 while (k > 0) {40 c = xr >= yr;41 q = (q << 1) | c;42 xr -= c ? yr : 0.0f;43 xr += xr;44 --k;45 }46 47 c = xr > yr;48 q = (q << 1) | c;49 xr -= c ? yr : 0.0f;50 51 int lt = ex < ey;52 53 q = lt ? 0 : q;54 xr = lt ? xa : xr;55 yr = lt ? ya : yr;56 57 c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & ((q & 0x1) == 0x1));58 xr -= c ? yr : 0.0f;59 q += c;60 61 float s = __clc_as_float(ey << EXPSHIFTBITS_SP32);62 xr *= lt ? 1.0f : s;63 64 c = ax == ay;65 xr = c ? 0.0f : xr;66 67 xr = __clc_as_float(sx ^ __clc_as_int(xr));68 69 c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 |70 ay == 0;71 xr = c ? __clc_as_float(QNANBITPATT_SP32) : xr;72 73 return xr;74}75 76#define __CLC_FLOAT_ONLY77#define __CLC_FUNCTION __clc_remainder78#define __CLC_BODY <clc/shared/binary_def_scalarize.inc>79#include <clc/math/gentype.inc>80#undef __CLC_FUNCTION81 82#ifdef cl_khr_fp6483 84#pragma OPENCL EXTENSION cl_khr_fp64 : enable85 86_CLC_DEF _CLC_OVERLOAD double __clc_remainder(double x, double y) {87 ulong ux = __clc_as_ulong(x);88 ulong ax = ux & ~SIGNBIT_DP64;89 ulong xsgn = ux ^ ax;90 double dx = __clc_as_double(ax);91 int xexp = __clc_convert_int(ax >> EXPSHIFTBITS_DP64);92 int xexp1 = 11 - (int)__clc_clz(ax & MANTBITS_DP64);93 xexp1 = xexp < 1 ? xexp1 : xexp;94 95 ulong uy = __clc_as_ulong(y);96 ulong ay = uy & ~SIGNBIT_DP64;97 double dy = __clc_as_double(ay);98 int yexp = __clc_convert_int(ay >> EXPSHIFTBITS_DP64);99 int yexp1 = 11 - (int)__clc_clz(ay & MANTBITS_DP64);100 yexp1 = yexp < 1 ? yexp1 : yexp;101 102 int qsgn = ((ux ^ uy) & SIGNBIT_DP64) == 0UL ? 1 : -1;103 104 // First assume |x| > |y|105 106 // Set ntimes to the number of times we need to do a107 // partial remainder. If the exponent of x is an exact multiple108 // of 53 larger than the exponent of y, and the mantissa of x is109 // less than the mantissa of y, ntimes will be one too large110 // but it doesn't matter - it just means that we'll go round111 // the loop below one extra time.112 int ntimes = __clc_max(0, (xexp1 - yexp1) / 53);113 double w = __clc_ldexp(dy, ntimes * 53);114 w = ntimes == 0 ? dy : w;115 double scale = ntimes == 0 ? 1.0 : 0x1.0p-53;116 117 // Each time round the loop we compute a partial remainder.118 // This is done by subtracting a large multiple of w119 // from x each time, where w is a scaled up version of y.120 // The subtraction must be performed exactly in quad121 // precision, though the result at each stage can122 // fit exactly in a double precision number.123 int i;124 double t, v, p, pp;125 126 for (i = 0; i < ntimes; i++) {127 // Compute integral multiplier128 t = __clc_trunc(dx / w);129 130 // Compute w * t in quad precision131 p = w * t;132 pp = __clc_fma(w, t, -p);133 134 // Subtract w * t from dx135 v = dx - p;136 dx = v + (((dx - v) - p) - pp);137 138 // If t was one too large, dx will be negative. Add back one w.139 dx += dx < 0.0 ? w : 0.0;140 141 // Scale w down by 2^(-53) for the next iteration142 w *= scale;143 }144 145 // One more time146 // Variable todd says whether the integer t is odd or not147 t = __clc_floor(dx / w);148 long lt = (long)t;149 int todd = lt & 1;150 151 p = w * t;152 pp = __clc_fma(w, t, -p);153 v = dx - p;154 dx = v + (((dx - v) - p) - pp);155 i = dx < 0.0;156 todd ^= i;157 dx += i ? w : 0.0;158 159 // At this point, dx lies in the range [0,dy)160 161 // For the fmod function, we're done apart from setting the correct sign.162 //163 // For the remainder function, we need to adjust dx164 // so that it lies in the range (-y/2, y/2] by carefully165 // subtracting w (== dy == y) if necessary. The rigmarole166 // with todd is to get the correct sign of the result167 // when x/y lies exactly half way between two integers,168 // when we need to choose the even integer.169 170 int al = (2.0 * dx > w) | (todd & (2.0 * dx == w));171 double dxl = dx - (al ? w : 0.0);172 173 int ag = (dx > 0.5 * w) | (todd & (dx == 0.5 * w));174 double dxg = dx - (ag ? w : 0.0);175 176 dx = dy < 0x1.0p+1022 ? dxl : dxg;177 178 double ret = __clc_as_double(xsgn ^ __clc_as_ulong(dx));179 dx = __clc_as_double(ax);180 181 // Now handle |x| == |y|182 int c = dx == dy;183 t = __clc_as_double(xsgn);184 ret = c ? t : ret;185 186 // Next, handle |x| < |y|187 c = dx < dy;188 ret = c ? x : ret;189 190 c &= (yexp<1023 & 2.0 * dx> dy) | (dx > 0.5 * dy);191 // we could use a conversion here instead since qsgn = +-1192 p = qsgn == 1 ? -1.0 : 1.0;193 t = __clc_fma(y, p, x);194 ret = c ? t : ret;195 196 // We don't need anything special for |x| == 0197 198 // |y| is 0199 c = dy == 0.0;200 ret = c ? __clc_as_double(QNANBITPATT_DP64) : ret;201 202 // y is +-Inf, NaN203 c = yexp > BIASEDEMAX_DP64;204 t = y == y ? x : y;205 ret = c ? t : ret;206 207 // x is +=Inf, NaN208 c = xexp > BIASEDEMAX_DP64;209 ret = c ? __clc_as_double(QNANBITPATT_DP64) : ret;210 211 return ret;212}213 214#define __CLC_DOUBLE_ONLY215#define __CLC_FUNCTION __clc_remainder216#define __CLC_BODY <clc/shared/binary_def_scalarize.inc>217#include <clc/math/gentype.inc>218#undef __CLC_FUNCTION219 220#endif221 222#ifdef cl_khr_fp16223 224#pragma OPENCL EXTENSION cl_khr_fp16 : enable225 226// Forward the half version of this builtin onto the float one227#define __CLC_HALF_ONLY228#define __CLC_FUNCTION __clc_remainder229#define __CLC_BODY <clc/math/binary_def_via_fp32.inc>230#include <clc/math/gentype.inc>231 232#endif233