libclc/generic/lib/math/clc_powr.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_fabs.h>
  26 #include <clc/math/clc_mad.h>
  27
  28 #include "config.h"
  29 #include "math.h"
  30 #include "tables.h"
  31
  32 // compute pow using log and exp
  33 // x^y = exp(y * log(x))
  34 //
  35 // we take care not to lose precision in the intermediate steps
  36 //
  37 // When computing log, calculate it in splits,
  38 //
  39 // r = f * (p_invead + p_inv_tail)
  40 // r = rh + rt
  41 //
  42 // calculate log polynomial using r, in end addition, do
  43 // poly = poly + ((rh-r) + rt)
  44 //
  45 // lth = -r
  46 // ltt = ((xexp * log2_t) - poly) + logT
  47 // lt = lth + ltt
  48 //
  49 // lh = (xexp * log2_h) + logH
  50 // l = lh + lt
  51 //
  52 // Calculate final log answer as gh and gt,
  53 // gh = l & higher-half bits
  54 // gt = (((ltt - (lt - lth)) + ((lh - l) + lt)) + (l - gh))
  55 //
  56 // yh = y & higher-half bits
  57 // yt = y - yh
  58 //
  59 // Before entering computation of exp,
  60 // vs = ((yt*gt + yt*gh) + yh*gt)
  61 // v = vs + yh*gh
  62 // vt = ((yh*gh - v) + vs)
  63 //
  64 // In calculation of exp, add vt to r that is used for poly
  65 // At the end of exp, do
  66 // ((((expT * poly) + expT) + expH*poly) + expH)
  67
  68 _CLC_DEF _CLC_OVERLOAD float __clc_powr(float x, float y) {
  69   int ix = as_int(x);
  70   int ax = ix & EXSIGNBIT_SP32;
  71   int xpos = ix == ax;
  72
  73   int iy = as_int(y);
  74   int ay = iy & EXSIGNBIT_SP32;
  75   int ypos = iy == ay;
  76
  77   // Extra precise log calculation
  78   // First handle case that x is close to 1
  79   float r = 1.0f - as_float(ax);
  80   int near1 = __clc_fabs(r) < 0x1.0p-4f;
  81   float r2 = r * r;
  82
  83   // Coefficients are just 1/3, 1/4, 1/5 and 1/6
  84   float poly = __clc_mad(
  85       r,
  86       __clc_mad(r,
  87                 __clc_mad(r, __clc_mad(r, 0x1.24924ap-3f, 0x1.555556p-3f),
  88                           0x1.99999ap-3f),
  89                 0x1.000000p-2f),
  90       0x1.555556p-2f);
  91
  92   poly *= r2 * r;
  93
  94   float lth_near1 = -r2 * 0.5f;
  95   float ltt_near1 = -poly;
  96   float lt_near1 = lth_near1 + ltt_near1;
  97   float lh_near1 = -r;
  98   float l_near1 = lh_near1 + lt_near1;
  99
 100   // Computations for x not near 1
 101   int m = (int)(ax >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
 102   float mf = (float)m;
 103   int ixs = as_int(as_float(ax | 0x3f800000) - 1.0f);
 104   float mfs = (float)((ixs >> EXPSHIFTBITS_SP32) - 253);
 105   int c = m == -127;
 106   int ixn = c ? ixs : ax;
 107   float mfn = c ? mfs : mf;
 108
 109   int indx = (ixn & 0x007f0000) + ((ixn & 0x00008000) << 1);
 110
 111   // F - Y
 112   float f = as_float(0x3f000000 | indx) -
 113             as_float(0x3f000000 | (ixn & MANTBITS_SP32));
 114
 115   indx = indx >> 16;
 116   float2 tv = USE_TABLE(log_inv_tbl_ep, indx);
 117   float rh = f * tv.s0;
 118   float rt = f * tv.s1;
 119   r = rh + rt;
 120
 121   poly = __clc_mad(r, __clc_mad(r, 0x1.0p-2f, 0x1.555556p-2f), 0x1.0p-1f) *
 122          (r * r);
 123   poly += (rh - r) + rt;
 124
 125   const float LOG2_HEAD = 0x1.62e000p-1f;  // 0.693115234
 126   const float LOG2_TAIL = 0x1.0bfbe8p-15f; // 0.0000319461833
 127   tv = USE_TABLE(loge_tbl, indx);
 128   float lth = -r;
 129   float ltt = __clc_mad(mfn, LOG2_TAIL, -poly) + tv.s1;
 130   float lt = lth + ltt;
 131   float lh = __clc_mad(mfn, LOG2_HEAD, tv.s0);
 132   float l = lh + lt;
 133
 134   // Select near 1 or not
 135   lth = near1 ? lth_near1 : lth;
 136   ltt = near1 ? ltt_near1 : ltt;
 137   lt = near1 ? lt_near1 : lt;
 138   lh = near1 ? lh_near1 : lh;
 139   l = near1 ? l_near1 : l;
 140
 141   float gh = as_float(as_int(l) & 0xfffff000);
 142   float gt = ((ltt - (lt - lth)) + ((lh - l) + lt)) + (l - gh);
 143
 144   float yh = as_float(iy & 0xfffff000);
 145
 146   float yt = y - yh;
 147
 148   float ylogx_s = __clc_mad(gt, yh, __clc_mad(gh, yt, yt * gt));
 149   float ylogx = __clc_mad(yh, gh, ylogx_s);
 150   float ylogx_t = __clc_mad(yh, gh, -ylogx) + ylogx_s;
 151
 152   // Extra precise exp of ylogx
 153   // 64/log2 : 92.332482616893657
 154   const float R_64_BY_LOG2 = 0x1.715476p+6f;
 155   int n = convert_int(ylogx * R_64_BY_LOG2);
 156   float nf = (float)n;
 157
 158   int j = n & 0x3f;
 159   m = n >> 6;
 160   int m2 = m << EXPSHIFTBITS_SP32;
 161   // log2/64 lead: 0.0108032227
 162   const float R_LOG2_BY_64_LD = 0x1.620000p-7f;
 163   // log2/64 tail: 0.0000272020388
 164   const float R_LOG2_BY_64_TL = 0x1.c85fdep-16f;
 165   r = __clc_mad(nf, -R_LOG2_BY_64_TL, __clc_mad(nf, -R_LOG2_BY_64_LD, ylogx)) +
 166       ylogx_t;
 167
 168   // Truncated Taylor series for e^r
 169   poly = __clc_mad(__clc_mad(__clc_mad(r, 0x1.555556p-5f, 0x1.555556p-3f), r,
 170                              0x1.000000p-1f),
 171                    r * r, r);
 172
 173   tv = USE_TABLE(exp_tbl_ep, j);
 174
 175   float expylogx =
 176       __clc_mad(tv.s0, poly, __clc_mad(tv.s1, poly, tv.s1)) + tv.s0;
 177   float sexpylogx = expylogx * as_float(0x1 << (m + 149));
 178   float texpylogx = as_float(as_int(expylogx) + m2);
 179   expylogx = m < -125 ? sexpylogx : texpylogx;
 180
 181   // Result is +-Inf if (ylogx + ylogx_t) > 128*log2
 182   expylogx = ((ylogx > 0x1.62e430p+6f) |
 183               (ylogx == 0x1.62e430p+6f & ylogx_t > -0x1.05c610p-22f))
 184                  ? as_float(PINFBITPATT_SP32)
 185                  : expylogx;
 186
 187   // Result is 0 if ylogx < -149*log2
 188   expylogx = ylogx < -0x1.9d1da0p+6f ? 0.0f : expylogx;
 189
 190   // Classify y:
 191   //   inty = 0 means not an integer.
 192   //   inty = 1 means odd integer.
 193   //   inty = 2 means even integer.
 194
 195   int yexp = (int)(ay >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32 + 1;
 196   int mask = (1 << (24 - yexp)) - 1;
 197   int yodd = ((iy >> (24 - yexp)) & 0x1) != 0;
 198   int inty = yodd ? 1 : 2;
 199   inty = (iy & mask) != 0 ? 0 : inty;
 200   inty = yexp < 1 ? 0 : inty;
 201   inty = yexp > 24 ? 2 : inty;
 202
 203   float signval = as_float((as_uint(expylogx) ^ SIGNBIT_SP32));
 204   expylogx = ((inty == 1) & !xpos) ? signval : expylogx;
 205   int ret = as_int(expylogx);
 206
 207   // Corner case handling
 208   ret = ax < 0x3f800000 & iy == NINFBITPATT_SP32 ? PINFBITPATT_SP32 : ret;
 209   ret = ax < 0x3f800000 & iy == PINFBITPATT_SP32 ? 0 : ret;
 210   ret = ax == 0x3f800000 & ay < PINFBITPATT_SP32 ? 0x3f800000 : ret;
 211   ret = ax == 0x3f800000 & ay == PINFBITPATT_SP32 ? QNANBITPATT_SP32 : ret;
 212   ret = ax > 0x3f800000 & iy == NINFBITPATT_SP32 ? 0 : ret;
 213   ret = ax > 0x3f800000 & iy == PINFBITPATT_SP32 ? PINFBITPATT_SP32 : ret;
 214   ret = ((ix < PINFBITPATT_SP32) & (ay == 0)) ? 0x3f800000 : ret;
 215   ret = ((ax == PINFBITPATT_SP32) & !ypos) ? 0 : ret;
 216   ret = ((ax == PINFBITPATT_SP32) & ypos) ? PINFBITPATT_SP32 : ret;
 217   ret = ((ax == PINFBITPATT_SP32) & (iy == PINFBITPATT_SP32)) ? PINFBITPATT_SP32
 218                                                               : ret;
 219   ret = ((ax == PINFBITPATT_SP32) & (ay == 0)) ? QNANBITPATT_SP32 : ret;
 220   ret = ((ax == 0) & !ypos) ? PINFBITPATT_SP32 : ret;
 221   ret = ((ax == 0) & ypos) ? 0 : ret;
 222   ret = ((ax == 0) & (ay == 0)) ? QNANBITPATT_SP32 : ret;
 223   ret = ((ax != 0) & !xpos) ? QNANBITPATT_SP32 : ret;
 224   ret = ax > PINFBITPATT_SP32 ? ix : ret;
 225   ret = ay > PINFBITPATT_SP32 ? iy : ret;
 226
 227   return as_float(ret);
 228 }
 229 _CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_powr, float, float)
 230
 231 #ifdef cl_khr_fp64
 232 _CLC_DEF _CLC_OVERLOAD double __clc_powr(double x, double y) {
 233   const double real_log2_tail = 5.76999904754328540596e-08;
 234   const double real_log2_lead = 6.93147122859954833984e-01;
 235
 236   long ux = as_long(x);
 237   long ax = ux & (~SIGNBIT_DP64);
 238   int xpos = ax == ux;
 239
 240   long uy = as_long(y);
 241   long ay = uy & (~SIGNBIT_DP64);
 242   int ypos = ay == uy;
 243
 244   // Extended precision log
 245   double v, vt;
 246   {
 247     int exp = (int)(ax >> 52) - 1023;
 248     int mask_exp_1023 = exp == -1023;
 249     double xexp = (double)exp;
 250     long mantissa = ax & 0x000FFFFFFFFFFFFFL;
 251
 252     long temp_ux = as_long(as_double(0x3ff0000000000000L | mantissa) - 1.0);
 253     exp = ((temp_ux & 0x7FF0000000000000L) >> 52) - 2045;
 254     double xexp1 = (double)exp;
 255     long mantissa1 = temp_ux & 0x000FFFFFFFFFFFFFL;
 256
 257     xexp = mask_exp_1023 ? xexp1 : xexp;
 258     mantissa = mask_exp_1023 ? mantissa1 : mantissa;
 259
 260     long rax = (mantissa & 0x000ff00000000000) +
 261                ((mantissa & 0x0000080000000000) << 1);
 262     int index = rax >> 44;
 263
 264     double F = as_double(rax | 0x3FE0000000000000L);
 265     double Y = as_double(mantissa | 0x3FE0000000000000L);
 266     double f = F - Y;
 267     double2 tv = USE_TABLE(log_f_inv_tbl, index);
 268     double log_h = tv.s0;
 269     double log_t = tv.s1;
 270     double f_inv = (log_h + log_t) * f;
 271     double r1 = as_double(as_long(f_inv) & 0xfffffffff8000000L);
 272     double r2 = fma(-F, r1, f) * (log_h + log_t);
 273     double r = r1 + r2;
 274
 275     double poly = fma(
 276         r, fma(r, fma(r, fma(r, 1.0 / 7.0, 1.0 / 6.0), 1.0 / 5.0), 1.0 / 4.0),
 277         1.0 / 3.0);
 278     poly = poly * r * r * r;
 279
 280     double hr1r1 = 0.5 * r1 * r1;
 281     double poly0h = r1 + hr1r1;
 282     double poly0t = r1 - poly0h + hr1r1;
 283     poly = fma(r1, r2, fma(0.5 * r2, r2, poly)) + r2 + poly0t;
 284
 285     tv = USE_TABLE(powlog_tbl, index);
 286     log_h = tv.s0;
 287     log_t = tv.s1;
 288
 289     double resT_t = fma(xexp, real_log2_tail, +log_t) - poly;
 290     double resT = resT_t - poly0h;
 291     double resH = fma(xexp, real_log2_lead, log_h);
 292     double resT_h = poly0h;
 293
 294     double H = resT + resH;
 295     double H_h = as_double(as_long(H) & 0xfffffffff8000000L);
 296     double T = (resH - H + resT) + (resT_t - (resT + resT_h)) + (H - H_h);
 297     H = H_h;
 298
 299     double y_head = as_double(uy & 0xfffffffff8000000L);
 300     double y_tail = y - y_head;
 301
 302     double temp = fma(y_tail, H, fma(y_head, T, y_tail * T));
 303     v = fma(y_head, H, temp);
 304     vt = fma(y_head, H, -v) + temp;
 305   }
 306
 307   // Now calculate exp of (v,vt)
 308
 309   double expv;
 310   {
 311     const double max_exp_arg = 709.782712893384;
 312     const double min_exp_arg = -745.1332191019411;
 313     const double sixtyfour_by_lnof2 = 92.33248261689366;
 314     const double lnof2_by_64_head = 0.010830424260348081;
 315     const double lnof2_by_64_tail = -4.359010638708991e-10;
 316
 317     double temp = v * sixtyfour_by_lnof2;
 318     int n = (int)temp;
 319     double dn = (double)n;
 320     int j = n & 0x0000003f;
 321     int m = n >> 6;
 322
 323     double2 tv = USE_TABLE(two_to_jby64_ep_tbl, j);
 324     double f1 = tv.s0;
 325     double f2 = tv.s1;
 326     double f = f1 + f2;
 327
 328     double r1 = fma(dn, -lnof2_by_64_head, v);
 329     double r2 = dn * lnof2_by_64_tail;
 330     double r = (r1 + r2) + vt;
 331
 332     double q = fma(
 333         r,
 334         fma(r,
 335             fma(r,
 336                 fma(r, 1.38889490863777199667e-03, 8.33336798434219616221e-03),
 337                 4.16666666662260795726e-02),
 338             1.66666666665260878863e-01),
 339         5.00000000000000008883e-01);
 340     q = fma(r * r, q, r);
 341
 342     expv = fma(f, q, f2) + f1;
 343     expv = ldexp(expv, m);
 344
 345     expv = v > max_exp_arg ? as_double(0x7FF0000000000000L) : expv;
 346     expv = v < min_exp_arg ? 0.0 : expv;
 347   }
 348
 349   // See whether y is an integer.
 350   // inty = 0 means not an integer.
 351   // inty = 1 means odd integer.
 352   // inty = 2 means even integer.
 353
 354   int inty;
 355   {
 356     int yexp = (int)(ay >> EXPSHIFTBITS_DP64) - EXPBIAS_DP64 + 1;
 357     inty = yexp < 1 ? 0 : 2;
 358     inty = yexp > 53 ? 2 : inty;
 359     long mask = (1L << (53 - yexp)) - 1L;
 360     int inty1 = (((ay & ~mask) >> (53 - yexp)) & 1L) == 1L ? 1 : 2;
 361     inty1 = (ay & mask) != 0 ? 0 : inty1;
 362     inty = !(yexp < 1) & !(yexp > 53) ? inty1 : inty;
 363   }
 364
 365   expv *= ((inty == 1) & !xpos) ? -1.0 : 1.0;
 366
 367   long ret = as_long(expv);
 368
 369   // Now all the edge cases
 370   ret = ax < 0x3ff0000000000000L & uy == NINFBITPATT_DP64 ? PINFBITPATT_DP64
 371                                                           : ret;
 372   ret = ax < 0x3ff0000000000000L & uy == PINFBITPATT_DP64 ? 0L : ret;
 373   ret = ax == 0x3ff0000000000000L & ay < PINFBITPATT_DP64 ? 0x3ff0000000000000L
 374                                                           : ret;
 375   ret = ax == 0x3ff0000000000000L & ay == PINFBITPATT_DP64 ? QNANBITPATT_DP64
 376                                                            : ret;
 377   ret = ax > 0x3ff0000000000000L & uy == NINFBITPATT_DP64 ? 0L : ret;
 378   ret = ax > 0x3ff0000000000000L & uy == PINFBITPATT_DP64 ? PINFBITPATT_DP64
 379                                                           : ret;
 380   ret = ux < PINFBITPATT_DP64 & ay == 0L ? 0x3ff0000000000000L : ret;
 381   ret = ((ax == PINFBITPATT_DP64) & !ypos) ? 0L : ret;
 382   ret = ((ax == PINFBITPATT_DP64) & ypos) ? PINFBITPATT_DP64 : ret;
 383   ret = ((ax == PINFBITPATT_DP64) & (uy == PINFBITPATT_DP64)) ? PINFBITPATT_DP64
 384                                                               : ret;
 385   ret = ((ax == PINFBITPATT_DP64) & (ay == 0L)) ? QNANBITPATT_DP64 : ret;
 386   ret = ((ax == 0L) & !ypos) ? PINFBITPATT_DP64 : ret;
 387   ret = ((ax == 0L) & ypos) ? 0L : ret;
 388   ret = ((ax == 0L) & (ay == 0L)) ? QNANBITPATT_DP64 : ret;
 389   ret = ((ax != 0L) & !xpos) ? QNANBITPATT_DP64 : ret;
 390   ret = ax > PINFBITPATT_DP64 ? ux : ret;
 391   ret = ay > PINFBITPATT_DP64 ? uy : ret;
 392
 393   return as_double(ret);
 394 }
 395 _CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, double, __clc_powr, double, double)
 396 #endif