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