1/* 2 * Copyright (c) 2014 Advanced Micro Devices, Inc. 3 * 4 * Permission is hereby granted, free of charge, to any person obtaining a copy 5 * of this software and associated documentation files (the "Software"), to deal 6 * in the Software without restriction, including without limitation the rights 7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 8 * copies of the Software, and to permit persons to whom the Software is 9 * furnished to do so, subject to the following conditions: 10 * 11 * The above copyright notice and this permission notice shall be included in 12 * all copies or substantial portions of the Software. 13 * 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 20 * THE SOFTWARE. 21 */ 22 23#include <clc/clc.h> 24 25#include <math/clc_remainder.h> 26#include "../clcmacro.h" 27#include "config.h" 28#include "math.h" 29 30_CLC_DEF _CLC_OVERLOAD float __clc_fmod(float x, float y) 31{ 32 int ux = as_int(x); 33 int ax = ux & EXSIGNBIT_SP32; 34 float xa = as_float(ax); 35 int sx = ux ^ ax; 36 int ex = ax >> EXPSHIFTBITS_SP32; 37 38 int uy = as_int(y); 39 int ay = uy & EXSIGNBIT_SP32; 40 float ya = as_float(ay); 41 int ey = ay >> EXPSHIFTBITS_SP32; 42 43 float xr = as_float(0x3f800000 | (ax & 0x007fffff)); 44 float yr = as_float(0x3f800000 | (ay & 0x007fffff)); 45 int c; 46 int k = ex - ey; 47 48 while (k > 0) { 49 c = xr >= yr; 50 xr -= c ? yr : 0.0f; 51 xr += xr; 52 --k; 53 } 54 55 c = xr >= yr; 56 xr -= c ? yr : 0.0f; 57 58 int lt = ex < ey; 59 60 xr = lt ? xa : xr; 61 yr = lt ? ya : yr; 62 63 64 float s = as_float(ey << EXPSHIFTBITS_SP32); 65 xr *= lt ? 1.0f : s; 66 67 c = ax == ay; 68 xr = c ? 0.0f : xr; 69 70 xr = as_float(sx ^ as_int(xr)); 71 72 c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | ay == 0; 73 xr = c ? as_float(QNANBITPATT_SP32) : xr; 74 75 return xr; 76 77} 78_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_fmod, float, float); 79 80#ifdef cl_khr_fp64 81_CLC_DEF _CLC_OVERLOAD double __clc_fmod(double x, double y) 82{ 83 ulong ux = as_ulong(x); 84 ulong ax = ux & ~SIGNBIT_DP64; 85 ulong xsgn = ux ^ ax; 86 double dx = as_double(ax); 87 int xexp = convert_int(ax >> EXPSHIFTBITS_DP64); 88 int xexp1 = 11 - (int) clz(ax & MANTBITS_DP64); 89 xexp1 = xexp < 1 ? xexp1 : xexp; 90 91 ulong uy = as_ulong(y); 92 ulong ay = uy & ~SIGNBIT_DP64; 93 double dy = as_double(ay); 94 int yexp = convert_int(ay >> EXPSHIFTBITS_DP64); 95 int yexp1 = 11 - (int) clz(ay & MANTBITS_DP64); 96 yexp1 = yexp < 1 ? yexp1 : yexp; 97 98 // First assume |x| > |y| 99 100 // Set ntimes to the number of times we need to do a 101 // partial remainder. If the exponent of x is an exact multiple 102 // of 53 larger than the exponent of y, and the mantissa of x is 103 // less than the mantissa of y, ntimes will be one too large 104 // but it doesn't matter - it just means that we'll go round 105 // the loop below one extra time. 106 int ntimes = max(0, (xexp1 - yexp1) / 53); 107 double w = ldexp(dy, ntimes * 53); 108 w = ntimes == 0 ? dy : w; 109 double scale = ntimes == 0 ? 1.0 : 0x1.0p-53; 110 111 // Each time round the loop we compute a partial remainder. 112 // This is done by subtracting a large multiple of w 113 // from x each time, where w is a scaled up version of y. 114 // The subtraction must be performed exactly in quad 115 // precision, though the result at each stage can 116 // fit exactly in a double precision number. 117 int i; 118 double t, v, p, pp; 119 120 for (i = 0; i < ntimes; i++) { 121 // Compute integral multiplier 122 t = trunc(dx / w); 123 124 // Compute w * t in quad precision 125 p = w * t; 126 pp = fma(w, t, -p); 127 128 // Subtract w * t from dx 129 v = dx - p; 130 dx = v + (((dx - v) - p) - pp); 131 132 // If t was one too large, dx will be negative. Add back one w. 133 dx += dx < 0.0 ? w : 0.0; 134 135 // Scale w down by 2^(-53) for the next iteration 136 w *= scale; 137 } 138 139 // One more time 140 // Variable todd says whether the integer t is odd or not 141 t = floor(dx / w); 142 long lt = (long)t; 143 int todd = lt & 1; 144 145 p = w * t; 146 pp = fma(w, t, -p); 147 v = dx - p; 148 dx = v + (((dx - v) - p) - pp); 149 i = dx < 0.0; 150 todd ^= i; 151 dx += i ? w : 0.0; 152 153 // At this point, dx lies in the range [0,dy) 154 double ret = as_double(xsgn ^ as_ulong(dx)); 155 dx = as_double(ax); 156 157 // Now handle |x| == |y| 158 int c = dx == dy; 159 t = as_double(xsgn); 160 ret = c ? t : ret; 161 162 // Next, handle |x| < |y| 163 c = dx < dy; 164 ret = c ? x : ret; 165 166 // We don't need anything special for |x| == 0 167 168 // |y| is 0 169 c = dy == 0.0; 170 ret = c ? as_double(QNANBITPATT_DP64) : ret; 171 172 // y is +-Inf, NaN 173 c = yexp > BIASEDEMAX_DP64; 174 t = y == y ? x : y; 175 ret = c ? t : ret; 176 177 // x is +=Inf, NaN 178 c = xexp > BIASEDEMAX_DP64; 179 ret = c ? as_double(QNANBITPATT_DP64) : ret; 180 181 return ret; 182} 183_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, double, __clc_fmod, double, double); 184#endif 185