197 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_trunc.h>16#include <clc/math/math.h>17#include <clc/shared/clc_max.h>18 19_CLC_DEF _CLC_OVERLOAD float __clc_fmod(float x, float y) {20 int ux = __clc_as_int(x);21 int ax = ux & EXSIGNBIT_SP32;22 float xa = __clc_as_float(ax);23 int sx = ux ^ ax;24 int ex = ax >> EXPSHIFTBITS_SP32;25 26 int uy = __clc_as_int(y);27 int ay = uy & EXSIGNBIT_SP32;28 float ya = __clc_as_float(ay);29 int ey = ay >> EXPSHIFTBITS_SP32;30 31 float xr = __clc_as_float(0x3f800000 | (ax & 0x007fffff));32 float yr = __clc_as_float(0x3f800000 | (ay & 0x007fffff));33 int c;34 int k = ex - ey;35 36 while (k > 0) {37 c = xr >= yr;38 xr -= c ? yr : 0.0f;39 xr += xr;40 --k;41 }42 43 c = xr >= yr;44 xr -= c ? yr : 0.0f;45 46 int lt = ex < ey;47 48 xr = lt ? xa : xr;49 yr = lt ? ya : yr;50 51 float s = __clc_as_float(ey << EXPSHIFTBITS_SP32);52 xr *= lt ? 1.0f : s;53 54 c = ax == ay;55 xr = c ? 0.0f : xr;56 57 xr = __clc_as_float(sx ^ __clc_as_int(xr));58 59 c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 |60 ay == 0;61 xr = c ? __clc_as_float(QNANBITPATT_SP32) : xr;62 63 return xr;64}65 66#define __CLC_FLOAT_ONLY67#define __CLC_FUNCTION __clc_fmod68#define __CLC_BODY <clc/shared/binary_def_scalarize.inc>69#include <clc/math/gentype.inc>70#undef __CLC_FUNCTION71 72#ifdef cl_khr_fp6473 74#pragma OPENCL EXTENSION cl_khr_fp64 : enable75 76_CLC_DEF _CLC_OVERLOAD double __clc_fmod(double x, double y) {77 ulong ux = __clc_as_ulong(x);78 ulong ax = ux & ~SIGNBIT_DP64;79 ulong xsgn = ux ^ ax;80 double dx = __clc_as_double(ax);81 int xexp = __clc_convert_int(ax >> EXPSHIFTBITS_DP64);82 int xexp1 = 11 - (int)__clc_clz(ax & MANTBITS_DP64);83 xexp1 = xexp < 1 ? xexp1 : xexp;84 85 ulong uy = __clc_as_ulong(y);86 ulong ay = uy & ~SIGNBIT_DP64;87 double dy = __clc_as_double(ay);88 int yexp = __clc_convert_int(ay >> EXPSHIFTBITS_DP64);89 int yexp1 = 11 - (int)__clc_clz(ay & MANTBITS_DP64);90 yexp1 = yexp < 1 ? yexp1 : yexp;91 92 // First assume |x| > |y|93 94 // Set ntimes to the number of times we need to do a95 // partial remainder. If the exponent of x is an exact multiple96 // of 53 larger than the exponent of y, and the mantissa of x is97 // less than the mantissa of y, ntimes will be one too large98 // but it doesn't matter - it just means that we'll go round99 // the loop below one extra time.100 int ntimes = __clc_max(0, (xexp1 - yexp1) / 53);101 double w = __clc_ldexp(dy, ntimes * 53);102 w = ntimes == 0 ? dy : w;103 double scale = ntimes == 0 ? 1.0 : 0x1.0p-53;104 105 // Each time round the loop we compute a partial remainder.106 // This is done by subtracting a large multiple of w107 // from x each time, where w is a scaled up version of y.108 // The subtraction must be performed exactly in quad109 // precision, though the result at each stage can110 // fit exactly in a double precision number.111 int i;112 double t, v, p, pp;113 114 for (i = 0; i < ntimes; i++) {115 // Compute integral multiplier116 t = __clc_trunc(dx / w);117 118 // Compute w * t in quad precision119 p = w * t;120 pp = __clc_fma(w, t, -p);121 122 // Subtract w * t from dx123 v = dx - p;124 dx = v + (((dx - v) - p) - pp);125 126 // If t was one too large, dx will be negative. Add back one w.127 dx += dx < 0.0 ? w : 0.0;128 129 // Scale w down by 2^(-53) for the next iteration130 w *= scale;131 }132 133 // One more time134 // Variable todd says whether the integer t is odd or not135 t = __clc_floor(dx / w);136 long lt = (long)t;137 int todd = lt & 1;138 139 p = w * t;140 pp = __clc_fma(w, t, -p);141 v = dx - p;142 dx = v + (((dx - v) - p) - pp);143 i = dx < 0.0;144 todd ^= i;145 dx += i ? w : 0.0;146 147 // At this point, dx lies in the range [0,dy)148 double ret = __clc_as_double(xsgn ^ __clc_as_ulong(dx));149 dx = __clc_as_double(ax);150 151 // Now handle |x| == |y|152 int c = dx == dy;153 t = __clc_as_double(xsgn);154 ret = c ? t : ret;155 156 // Next, handle |x| < |y|157 c = dx < dy;158 ret = c ? x : ret;159 160 // We don't need anything special for |x| == 0161 162 // |y| is 0163 c = dy == 0.0;164 ret = c ? __clc_as_double(QNANBITPATT_DP64) : ret;165 166 // y is +-Inf, NaN167 c = yexp > BIASEDEMAX_DP64;168 t = y == y ? x : y;169 ret = c ? t : ret;170 171 // x is +=Inf, NaN172 c = xexp > BIASEDEMAX_DP64;173 ret = c ? __clc_as_double(QNANBITPATT_DP64) : ret;174 175 return ret;176}177 178#define __CLC_DOUBLE_ONLY179#define __CLC_FUNCTION __clc_fmod180#define __CLC_BODY <clc/shared/binary_def_scalarize.inc>181#include <clc/math/gentype.inc>182#undef __CLC_FUNCTION183 184#endif185 186#ifdef cl_khr_fp16187 188#pragma OPENCL EXTENSION cl_khr_fp16 : enable189 190// Forward the half version of this builtin onto the float one191#define __CLC_HALF_ONLY192#define __CLC_FUNCTION __clc_fmod193#define __CLC_BODY <clc/math/binary_def_via_fp32.inc>194#include <clc/math/gentype.inc>195 196#endif197