libclc/generic/lib/math/clc_remquo.cl

   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/math/clc_floor.h>
  26 #include <clc/math/clc_trunc.h>
  27 #include <clc/shared/clc_max.h>
  28 #include <math/clc_remainder.h>
  29
  30 #include "config.h"
  31 #include "math.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)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)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