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_remainder(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 uint q = 0; 49 50 while (k > 0) { 51 c = xr >= yr; 52 q = (q << 1) | c; 53 xr -= c ? yr : 0.0f; 54 xr += xr; 55 --k; 56 } 57 58 c = xr > yr; 59 q = (q << 1) | c; 60 xr -= c ? yr : 0.0f; 61 62 int lt = ex < ey; 63 64 q = lt ? 0 : q; 65 xr = lt ? xa : xr; 66 yr = lt ? ya : yr; 67 68 c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & ((q & 0x1) == 0x1)); 69 xr -= c ? yr : 0.0f; 70 q += c; 71 72 float s = as_float(ey << EXPSHIFTBITS_SP32); 73 xr *= lt ? 1.0f : s; 74 75 c = ax == ay; 76 xr = c ? 0.0f : xr; 77 78 xr = as_float(sx ^ as_int(xr)); 79 80 c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | ay == 0; 81 xr = c ? as_float(QNANBITPATT_SP32) : xr; 82 83 return xr; 84 85} 86_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_remainder, float, float); 87 88#ifdef cl_khr_fp64 89_CLC_DEF _CLC_OVERLOAD double __clc_remainder(double x, double y) 90{ 91 ulong ux = as_ulong(x); 92 ulong ax = ux & ~SIGNBIT_DP64; 93 ulong xsgn = ux ^ ax; 94 double dx = as_double(ax); 95 int xexp = convert_int(ax >> EXPSHIFTBITS_DP64); 96 int xexp1 = 11 - (int) clz(ax & MANTBITS_DP64); 97 xexp1 = xexp < 1 ? xexp1 : xexp; 98 99 ulong uy = as_ulong(y); 100 ulong ay = uy & ~SIGNBIT_DP64; 101 double dy = as_double(ay); 102 int yexp = convert_int(ay >> EXPSHIFTBITS_DP64); 103 int yexp1 = 11 - (int) clz(ay & MANTBITS_DP64); 104 yexp1 = yexp < 1 ? yexp1 : yexp; 105 106 int qsgn = ((ux ^ uy) & SIGNBIT_DP64) == 0UL ? 1 : -1; 107 108 // First assume |x| > |y| 109 110 // Set ntimes to the number of times we need to do a 111 // partial remainder. If the exponent of x is an exact multiple 112 // of 53 larger than the exponent of y, and the mantissa of x is 113 // less than the mantissa of y, ntimes will be one too large 114 // but it doesn't matter - it just means that we'll go round 115 // the loop below one extra time. 116 int ntimes = max(0, (xexp1 - yexp1) / 53); 117 double w = ldexp(dy, ntimes * 53); 118 w = ntimes == 0 ? dy : w; 119 double scale = ntimes == 0 ? 1.0 : 0x1.0p-53; 120 121 // Each time round the loop we compute a partial remainder. 122 // This is done by subtracting a large multiple of w 123 // from x each time, where w is a scaled up version of y. 124 // The subtraction must be performed exactly in quad 125 // precision, though the result at each stage can 126 // fit exactly in a double precision number. 127 int i; 128 double t, v, p, pp; 129 130 for (i = 0; i < ntimes; i++) { 131 // Compute integral multiplier 132 t = trunc(dx / w); 133 134 // Compute w * t in quad precision 135 p = w * t; 136 pp = fma(w, t, -p); 137 138 // Subtract w * t from dx 139 v = dx - p; 140 dx = v + (((dx - v) - p) - pp); 141 142 // If t was one too large, dx will be negative. Add back one w. 143 dx += dx < 0.0 ? w : 0.0; 144 145 // Scale w down by 2^(-53) for the next iteration 146 w *= scale; 147 } 148 149 // One more time 150 // Variable todd says whether the integer t is odd or not 151 t = floor(dx / w); 152 long lt = (long)t; 153 int todd = lt & 1; 154 155 p = w * t; 156 pp = fma(w, t, -p); 157 v = dx - p; 158 dx = v + (((dx - v) - p) - pp); 159 i = dx < 0.0; 160 todd ^= i; 161 dx += i ? w : 0.0; 162 163 // At this point, dx lies in the range [0,dy) 164 165 // For the fmod function, we're done apart from setting the correct sign. 166 // 167 // For the remainder function, we need to adjust dx 168 // so that it lies in the range (-y/2, y/2] by carefully 169 // subtracting w (== dy == y) if necessary. The rigmarole 170 // with todd is to get the correct sign of the result 171 // when x/y lies exactly half way between two integers, 172 // when we need to choose the even integer. 173 174 int al = (2.0*dx > w) | (todd & (2.0*dx == w)); 175 double dxl = dx - (al ? w : 0.0); 176 177 int ag = (dx > 0.5*w) | (todd & (dx == 0.5*w)); 178 double dxg = dx - (ag ? w : 0.0); 179 180 dx = dy < 0x1.0p+1022 ? dxl : dxg; 181 182 double ret = as_double(xsgn ^ as_ulong(dx)); 183 dx = as_double(ax); 184 185 // Now handle |x| == |y| 186 int c = dx == dy; 187 t = as_double(xsgn); 188 ret = c ? t : ret; 189 190 // Next, handle |x| < |y| 191 c = dx < dy; 192 ret = c ? x : ret; 193 194 c &= (yexp < 1023 & 2.0*dx > dy) | (dx > 0.5*dy); 195 // we could use a conversion here instead since qsgn = +-1 196 p = qsgn == 1 ? -1.0 : 1.0; 197 t = fma(y, p, x); 198 ret = c ? t : ret; 199 200 // We don't need anything special for |x| == 0 201 202 // |y| is 0 203 c = dy == 0.0; 204 ret = c ? as_double(QNANBITPATT_DP64) : ret; 205 206 // y is +-Inf, NaN 207 c = yexp > BIASEDEMAX_DP64; 208 t = y == y ? x : y; 209 ret = c ? t : ret; 210 211 // x is +=Inf, NaN 212 c = xexp > BIASEDEMAX_DP64; 213 ret = c ? as_double(QNANBITPATT_DP64) : ret; 214 215 return ret; 216} 217_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, double, __clc_remainder, double, double); 218#endif 219