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