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#include <clc/clcmacro.h> 25#include <clc/integer/clc_clz.h> 26#include <clc/math/clc_floor.h> 27#include <clc/math/clc_subnormal_config.h> 28#include <clc/math/clc_trunc.h> 29#include <clc/math/math.h> 30#include <clc/shared/clc_max.h> 31#include <math/clc_remainder.h> 32 33_CLC_DEF _CLC_OVERLOAD float __clc_remquo(float x, float y, 34 __private int *quo) { 35 x = __clc_flush_denormal_if_not_supported(x); 36 y = __clc_flush_denormal_if_not_supported(y); 37 int ux = as_int(x); 38 int ax = ux & EXSIGNBIT_SP32; 39 float xa = as_float(ax); 40 int sx = ux ^ ax; 41 int ex = ax >> EXPSHIFTBITS_SP32; 42 43 int uy = as_int(y); 44 int ay = uy & EXSIGNBIT_SP32; 45 float ya = as_float(ay); 46 int sy = uy ^ ay; 47 int ey = ay >> EXPSHIFTBITS_SP32; 48 49 float xr = as_float(0x3f800000 | (ax & 0x007fffff)); 50 float yr = as_float(0x3f800000 | (ay & 0x007fffff)); 51 int c; 52 int k = ex - ey; 53 54 uint q = 0; 55 56 while (k > 0) { 57 c = xr >= yr; 58 q = (q << 1) | c; 59 xr -= c ? yr : 0.0f; 60 xr += xr; 61 --k; 62 } 63 64 c = xr > yr; 65 q = (q << 1) | c; 66 xr -= c ? yr : 0.0f; 67 68 int lt = ex < ey; 69 70 q = lt ? 0 : q; 71 xr = lt ? xa : xr; 72 yr = lt ? ya : yr; 73 74 c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & ((q & 0x1) == 0x1)); 75 xr -= c ? yr : 0.0f; 76 q += c; 77 78 float s = as_float(ey << EXPSHIFTBITS_SP32); 79 xr *= lt ? 1.0f : s; 80 81 int qsgn = sx == sy ? 1 : -1; 82 int quot = (q & 0x7f) * qsgn; 83 84 c = ax == ay; 85 quot = c ? qsgn : quot; 86 xr = c ? 0.0f : xr; 87 88 xr = as_float(sx ^ as_int(xr)); 89 90 c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | 91 ay == 0; 92 quot = c ? 0 : quot; 93 xr = c ? as_float(QNANBITPATT_SP32) : xr; 94 95 *quo = quot; 96 97 return xr; 98} 99// remquo signature is special, we don't have macro for this 100#define __VEC_REMQUO(TYPE, VEC_SIZE, HALF_VEC_SIZE) \ 101 _CLC_DEF _CLC_OVERLOAD TYPE##VEC_SIZE __clc_remquo( \ 102 TYPE##VEC_SIZE x, TYPE##VEC_SIZE y, __private int##VEC_SIZE *quo) { \ 103 int##HALF_VEC_SIZE lo, hi; \ 104 TYPE##VEC_SIZE ret; \ 105 ret.lo = __clc_remquo(x.lo, y.lo, &lo); \ 106 ret.hi = __clc_remquo(x.hi, y.hi, &hi); \ 107 (*quo).lo = lo; \ 108 (*quo).hi = hi; \ 109 return ret; \ 110 } 111 112#define __VEC3_REMQUO(TYPE) \ 113 _CLC_DEF _CLC_OVERLOAD TYPE##3 __clc_remquo(TYPE##3 x, TYPE##3 y, \ 114 __private int##3 * quo) { \ 115 int2 lo; \ 116 int hi; \ 117 TYPE##3 ret; \ 118 ret.s01 = __clc_remquo(x.s01, y.s01, &lo); \ 119 ret.s2 = __clc_remquo(x.s2, y.s2, &hi); \ 120 (*quo).s01 = lo; \ 121 (*quo).s2 = hi; \ 122 return ret; \ 123 } 124__VEC_REMQUO(float, 2, ) 125__VEC3_REMQUO(float) 126__VEC_REMQUO(float, 4, 2) 127__VEC_REMQUO(float, 8, 4) 128__VEC_REMQUO(float, 16, 8) 129 130#ifdef cl_khr_fp64 131_CLC_DEF _CLC_OVERLOAD double __clc_remquo(double x, double y, 132 __private int *pquo) { 133 ulong ux = as_ulong(x); 134 ulong ax = ux & ~SIGNBIT_DP64; 135 ulong xsgn = ux ^ ax; 136 double dx = as_double(ax); 137 int xexp = convert_int(ax >> EXPSHIFTBITS_DP64); 138 int xexp1 = 11 - (int)__clc_clz(ax & MANTBITS_DP64); 139 xexp1 = xexp < 1 ? xexp1 : xexp; 140 141 ulong uy = as_ulong(y); 142 ulong ay = uy & ~SIGNBIT_DP64; 143 double dy = as_double(ay); 144 int yexp = convert_int(ay >> EXPSHIFTBITS_DP64); 145 int yexp1 = 11 - (int)__clc_clz(ay & MANTBITS_DP64); 146 yexp1 = yexp < 1 ? yexp1 : yexp; 147 148 int qsgn = ((ux ^ uy) & SIGNBIT_DP64) == 0UL ? 1 : -1; 149 150 // First assume |x| > |y| 151 152 // Set ntimes to the number of times we need to do a 153 // partial remainder. If the exponent of x is an exact multiple 154 // of 53 larger than the exponent of y, and the mantissa of x is 155 // less than the mantissa of y, ntimes will be one too large 156 // but it doesn't matter - it just means that we'll go round 157 // the loop below one extra time. 158 int ntimes = __clc_max(0, (xexp1 - yexp1) / 53); 159 double w = ldexp(dy, ntimes * 53); 160 w = ntimes == 0 ? dy : w; 161 double scale = ntimes == 0 ? 1.0 : 0x1.0p-53; 162 163 // Each time round the loop we compute a partial remainder. 164 // This is done by subtracting a large multiple of w 165 // from x each time, where w is a scaled up version of y. 166 // The subtraction must be performed exactly in quad 167 // precision, though the result at each stage can 168 // fit exactly in a double precision number. 169 int i; 170 double t, v, p, pp; 171 172 for (i = 0; i < ntimes; i++) { 173 // Compute integral multiplier 174 t = __clc_trunc(dx / w); 175 176 // Compute w * t in quad precision 177 p = w * t; 178 pp = fma(w, t, -p); 179 180 // Subtract w * t from dx 181 v = dx - p; 182 dx = v + (((dx - v) - p) - pp); 183 184 // If t was one too large, dx will be negative. Add back one w. 185 dx += dx < 0.0 ? w : 0.0; 186 187 // Scale w down by 2^(-53) for the next iteration 188 w *= scale; 189 } 190 191 // One more time 192 // Variable todd says whether the integer t is odd or not 193 t = __clc_floor(dx / w); 194 long lt = (long)t; 195 int todd = lt & 1; 196 197 p = w * t; 198 pp = fma(w, t, -p); 199 v = dx - p; 200 dx = v + (((dx - v) - p) - pp); 201 i = dx < 0.0; 202 todd ^= i; 203 dx += i ? w : 0.0; 204 205 lt -= i; 206 207 // At this point, dx lies in the range [0,dy) 208 209 // For the remainder function, we need to adjust dx 210 // so that it lies in the range (-y/2, y/2] by carefully 211 // subtracting w (== dy == y) if necessary. The rigmarole 212 // with todd is to get the correct sign of the result 213 // when x/y lies exactly half way between two integers, 214 // when we need to choose the even integer. 215 216 int al = (2.0 * dx > w) | (todd & (2.0 * dx == w)); 217 double dxl = dx - (al ? w : 0.0); 218 219 int ag = (dx > 0.5 * w) | (todd & (dx == 0.5 * w)); 220 double dxg = dx - (ag ? w : 0.0); 221 222 dx = dy < 0x1.0p+1022 ? dxl : dxg; 223 lt += dy < 0x1.0p+1022 ? al : ag; 224 int quo = ((int)lt & 0x7f) * qsgn; 225 226 double ret = as_double(xsgn ^ as_ulong(dx)); 227 dx = as_double(ax); 228 229 // Now handle |x| == |y| 230 int c = dx == dy; 231 t = as_double(xsgn); 232 quo = c ? qsgn : quo; 233 ret = c ? t : ret; 234 235 // Next, handle |x| < |y| 236 c = dx < dy; 237 quo = c ? 0 : quo; 238 ret = c ? x : ret; 239 240 c &= (yexp<1023 & 2.0 * dx> dy) | (dx > 0.5 * dy); 241 quo = c ? qsgn : quo; 242 // we could use a conversion here instead since qsgn = +-1 243 p = qsgn == 1 ? -1.0 : 1.0; 244 t = fma(y, p, x); 245 ret = c ? t : ret; 246 247 // We don't need anything special for |x| == 0 248 249 // |y| is 0 250 c = dy == 0.0; 251 quo = c ? 0 : quo; 252 ret = c ? as_double(QNANBITPATT_DP64) : ret; 253 254 // y is +-Inf, NaN 255 c = yexp > BIASEDEMAX_DP64; 256 quo = c ? 0 : quo; 257 t = y == y ? x : y; 258 ret = c ? t : ret; 259 260 // x is +=Inf, NaN 261 c = xexp > BIASEDEMAX_DP64; 262 quo = c ? 0 : quo; 263 ret = c ? as_double(QNANBITPATT_DP64) : ret; 264 265 *pquo = quo; 266 return ret; 267} 268__VEC_REMQUO(double, 2, ) 269__VEC3_REMQUO(double) 270__VEC_REMQUO(double, 4, 2) 271__VEC_REMQUO(double, 8, 4) 272__VEC_REMQUO(double, 16, 8) 273#endif 274 275#ifdef cl_khr_fp16 276 277#pragma OPENCL EXTENSION cl_khr_fp16 : enable 278 279_CLC_OVERLOAD _CLC_DEF half __clc_remquo(half x, half y, __private int *pquo) { 280 return (half)__clc_remquo((float)x, (float)y, pquo); 281} 282__VEC_REMQUO(half, 2, ) 283__VEC3_REMQUO(half) 284__VEC_REMQUO(half, 4, 2) 285__VEC_REMQUO(half, 8, 4) 286__VEC_REMQUO(half, 16, 8) 287 288#endif 289