403 lines · plain
1//===----------------------------------------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8//9// Computes pow using log and exp10//11// x^y = exp(y * log(x))12//13// We take care not to lose precision in the intermediate steps.14//15// When computing log, calculate it in splits:16//17// r = f * (p_invead + p_inv_tail)18// r = rh + rt19//20// Calculate log polynomial using r, in end addition, do:21//22// poly = poly + ((rh-r) + rt)23//24// lth = -r25// ltt = ((xexp * log2_t) - poly) + logT26// lt = lth + ltt27//28// lh = (xexp * log2_h) + logH29// l = lh + lt30//31// Calculate final log answer as gh and gt:32//33// gh = l & higher-half bits34// gt = (((ltt - (lt - lth)) + ((lh - l) + lt)) + (l - gh))35//36// yh = y & higher-half bits37// yt = y - yh38//39// Before entering computation of exp:40//41// vs = ((yt*gt + yt*gh) + yh*gt)42// v = vs + yh*gh43// vt = ((yh*gh - v) + vs)44//45// In calculation of exp, add vt to r that is used for poly.46//47// At the end of exp, do:48//49// ((((expT * poly) + expT) + expH*poly) + expH)50//51//===----------------------------------------------------------------------===//52 53#if __CLC_FPSIZE == 3254 55_CLC_DEF _CLC_OVERLOAD __CLC_GENTYPE __clc_pown(__CLC_GENTYPE x,56 __CLC_INTN ny) {57 __CLC_GENTYPE y = __CLC_CONVERT_GENTYPE(ny);58 59 __CLC_GENTYPE absx = __clc_fabs(x);60 __CLC_INTN ix = __CLC_AS_INTN(x);61 __CLC_INTN ax = __CLC_AS_INTN(absx);62 __CLC_INTN xpos = ix == ax;63 64 __CLC_INTN iy = __CLC_AS_INTN(y);65 __CLC_INTN ay = __CLC_AS_INTN(__clc_fabs(y));66 __CLC_INTN ypos = iy == ay;67 68 // Extra precise log calculation69 // First handle case that x is close to 170 __CLC_GENTYPE r = 1.0f - absx;71 __CLC_INTN near1 = __clc_fabs(r) < 0x1.0p-4f;72 __CLC_GENTYPE r2 = r * r;73 74 // Coefficients are just 1/3, 1/4, 1/5 and 1/675 __CLC_GENTYPE poly = __clc_mad(76 r,77 __clc_mad(r,78 __clc_mad(r, __clc_mad(r, 0x1.24924ap-3f, 0x1.555556p-3f),79 0x1.99999ap-3f),80 0x1.000000p-2f),81 0x1.555556p-2f);82 83 poly *= r2 * r;84 85 __CLC_GENTYPE lth_near1 = -r2 * 0.5f;86 __CLC_GENTYPE ltt_near1 = -poly;87 __CLC_GENTYPE lt_near1 = lth_near1 + ltt_near1;88 __CLC_GENTYPE lh_near1 = -r;89 __CLC_GENTYPE l_near1 = lh_near1 + lt_near1;90 91 // Computations for x not near 192 __CLC_INTN m = __CLC_CONVERT_INTN(ax >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;93 __CLC_GENTYPE mf = __CLC_CONVERT_GENTYPE(m);94 __CLC_INTN ixs = __CLC_AS_INTN(__CLC_AS_GENTYPE(ax | 0x3f800000) - 1.0f);95 __CLC_GENTYPE mfs = __CLC_CONVERT_GENTYPE((ixs >> EXPSHIFTBITS_SP32) - 253);96 __CLC_INTN c = m == -127;97 __CLC_INTN ixn = c ? ixs : ax;98 __CLC_GENTYPE mfn = c ? mfs : mf;99 100 __CLC_INTN indx = (ixn & 0x007f0000) + ((ixn & 0x00008000) << 1);101 102 // F - Y103 __CLC_GENTYPE f = __CLC_AS_GENTYPE(0x3f000000 | indx) -104 __CLC_AS_GENTYPE(0x3f000000 | (ixn & MANTBITS_SP32));105 106 indx = indx >> 16;107 __CLC_GENTYPE rh = f * __CLC_USE_TABLE(log_inv_tbl_ep_head, indx);108 __CLC_GENTYPE rt = f * __CLC_USE_TABLE(log_inv_tbl_ep_tail, indx);109 r = rh + rt;110 111 poly = __clc_mad(r, __clc_mad(r, 0x1.0p-2f, 0x1.555556p-2f), 0x1.0p-1f) *112 (r * r);113 poly += (rh - r) + rt;114 115 const __CLC_GENTYPE LOG2_HEAD = 0x1.62e000p-1f; // 0.693115234116 const __CLC_GENTYPE LOG2_TAIL = 0x1.0bfbe8p-15f; // 0.0000319461833117 __CLC_GENTYPE logel = __CLC_USE_TABLE(loge_tbl_lo, indx);118 __CLC_GENTYPE logeh = __CLC_USE_TABLE(loge_tbl_hi, indx);119 __CLC_GENTYPE lth = -r;120 __CLC_GENTYPE ltt = __clc_mad(mfn, LOG2_TAIL, -poly) + logeh;121 __CLC_GENTYPE lt = lth + ltt;122 __CLC_GENTYPE lh = __clc_mad(mfn, LOG2_HEAD, logel);123 __CLC_GENTYPE l = lh + lt;124 125 // Select near 1 or not126 lth = near1 ? lth_near1 : lth;127 ltt = near1 ? ltt_near1 : ltt;128 lt = near1 ? lt_near1 : lt;129 lh = near1 ? lh_near1 : lh;130 l = near1 ? l_near1 : l;131 132 __CLC_GENTYPE gh = __CLC_AS_GENTYPE(__CLC_AS_UINTN(l) & 0xfffff000);133 __CLC_GENTYPE gt = ((ltt - (lt - lth)) + ((lh - l) + lt)) + (l - gh);134 135 __CLC_GENTYPE yh = __CLC_AS_GENTYPE(__CLC_AS_UINTN(iy) & 0xfffff000);136 137 __CLC_GENTYPE yt = __CLC_CONVERT_GENTYPE(ny - __CLC_CONVERT_INTN(yh));138 139 __CLC_GENTYPE ylogx_s = __clc_mad(gt, yh, __clc_mad(gh, yt, yt * gt));140 __CLC_GENTYPE ylogx = __clc_mad(yh, gh, ylogx_s);141 __CLC_GENTYPE ylogx_t = __clc_mad(yh, gh, -ylogx) + ylogx_s;142 143 // Extra precise exp of ylogx144 // 64/log2 : 92.332482616893657145 const __CLC_GENTYPE R_64_BY_LOG2 = 0x1.715476p+6f;146 __CLC_INTN n = __CLC_CONVERT_INTN(ylogx * R_64_BY_LOG2);147 __CLC_GENTYPE nf = __CLC_CONVERT_GENTYPE(n);148 149 __CLC_INTN j = n & 0x3f;150 m = n >> 6;151 __CLC_INTN m2 = m << EXPSHIFTBITS_SP32;152 153 // log2/64 lead: 0.0108032227154 const __CLC_GENTYPE R_LOG2_BY_64_LD = 0x1.620000p-7f;155 // log2/64 tail: 0.0000272020388156 const __CLC_GENTYPE R_LOG2_BY_64_TL = 0x1.c85fdep-16f;157 r = __clc_mad(nf, -R_LOG2_BY_64_TL, __clc_mad(nf, -R_LOG2_BY_64_LD, ylogx)) +158 ylogx_t;159 160 // Truncated Taylor series for e^r161 poly = __clc_mad(__clc_mad(__clc_mad(r, 0x1.555556p-5f, 0x1.555556p-3f), r,162 0x1.000000p-1f),163 r * r, r);164 165 __CLC_GENTYPE exp_head = __CLC_USE_TABLE(exp_tbl_ep_head, j);166 __CLC_GENTYPE exp_tail = __CLC_USE_TABLE(exp_tbl_ep_tail, j);167 168 __CLC_GENTYPE expylogx =169 __clc_mad(exp_head, poly, __clc_mad(exp_tail, poly, exp_tail)) + exp_head;170 __CLC_GENTYPE sexpylogx =171 expylogx * __CLC_AS_GENTYPE((__CLC_INTN)0x1 << (m + 149));172 __CLC_GENTYPE texpylogx = __CLC_AS_GENTYPE(__CLC_AS_INTN(expylogx) + m2);173 expylogx = m < -125 ? sexpylogx : texpylogx;174 175 // Result is +-Inf if (ylogx + ylogx_t) > 128*log2176 expylogx =177 __clc_select(expylogx, __CLC_AS_GENTYPE((__CLC_UINTN)PINFBITPATT_SP32),178 ylogx > 0x1.62e430p+6f ||179 (ylogx == 0x1.62e430p+6f && ylogx_t > -0x1.05c610p-22f));180 181 // Result is 0 if ylogx < -149*log2182 expylogx = ylogx < -0x1.9d1da0p+6f ? 0.0f : expylogx;183 184 // Classify y:185 // inty = 0 means not an integer.186 // inty = 1 means odd integer.187 // inty = 2 means even integer.188 189 __CLC_INTN inty = 2 - (ny & 1);190 191 __CLC_GENTYPE signval =192 __CLC_AS_GENTYPE((__CLC_AS_UINTN(expylogx) ^ SIGNBIT_SP32));193 expylogx = ((inty == 1) && !xpos) ? signval : expylogx;194 __CLC_INTN ret = __CLC_AS_INTN(expylogx);195 196 // Corner case handling197 __CLC_BIT_INTN x_is_ninf = ix == (__CLC_INTN)NINFBITPATT_SP32;198 199 __CLC_INTN xinf =200 xpos ? (__CLC_INTN)PINFBITPATT_SP32 : (__CLC_INTN)NINFBITPATT_SP32;201 ret = ((ax == 0) && !ypos && (inty == 1)) ? xinf : ret;202 ret = ((ax == 0) && !ypos && (inty == 2)) ? PINFBITPATT_SP32 : ret;203 ret = ((ax == 0) && ypos && (inty == 2)) ? 0 : ret;204 __CLC_INTN xzero = !xpos ? (__CLC_INTN)0x80000000 : (__CLC_INTN)0;205 ret = ((ax == 0) && ypos && (inty == 1)) ? xzero : ret;206 ret = (x_is_ninf && !ypos && (inty == 1)) ? (__CLC_INTN)0x80000000 : ret;207 ret = (x_is_ninf && !ypos && (inty != 1)) ? 0 : ret;208 ret = (x_is_ninf && ypos && (inty == 1)) ? (__CLC_INTN)NINFBITPATT_SP32 : ret;209 ret = (x_is_ninf && ypos && (inty != 1)) ? (__CLC_INTN)PINFBITPATT_SP32 : ret;210 ret = ((ix == PINFBITPATT_SP32) && !ypos) ? 0 : ret;211 ret = ((ix == PINFBITPATT_SP32) && ypos) ? (__CLC_INTN)PINFBITPATT_SP32 : ret;212 ret = ax > PINFBITPATT_SP32 ? ix : ret;213 ret = ny == 0 ? 0x3f800000 : ret;214 215 return __CLC_AS_GENTYPE(ret);216}217 218#elif __CLC_FPSIZE == 64219 220_CLC_DEF _CLC_OVERLOAD __CLC_GENTYPE __clc_pown(__CLC_GENTYPE x,221 __CLC_INTN ny) {222 const __CLC_GENTYPE real_log2_tail = 5.76999904754328540596e-08;223 const __CLC_GENTYPE real_log2_lead = 6.93147122859954833984e-01;224 225 __CLC_GENTYPE y = __CLC_CONVERT_GENTYPE(ny);226 227 __CLC_LONGN ux = __CLC_AS_LONGN(x);228 __CLC_LONGN ax = __CLC_AS_LONGN(__clc_fabs(x));229 __CLC_BIT_INTN xpos = ax == ux;230 231 __CLC_LONGN uy = __CLC_AS_LONGN(y);232 __CLC_LONGN ay = __CLC_AS_LONGN(__clc_fabs(y));233 __CLC_BIT_INTN ypos = ay == uy;234 235 // Extended precision log236 __CLC_GENTYPE v, vt;237 {238 __CLC_INTN exp = __CLC_CONVERT_INTN(ax >> 52) - 1023;239 __CLC_INTN mask_exp_1023 = exp == -1023;240 __CLC_GENTYPE xexp = __CLC_CONVERT_GENTYPE(exp);241 __CLC_LONGN mantissa = ax & 0x000FFFFFFFFFFFFFL;242 243 __CLC_LONGN temp_ux =244 __CLC_AS_LONGN(__CLC_AS_GENTYPE(0x3ff0000000000000L | mantissa) - 1.0);245 exp = __CLC_CONVERT_INTN((temp_ux & 0x7FF0000000000000L) >> 52) - 2045;246 __CLC_GENTYPE xexp1 = __CLC_CONVERT_GENTYPE(exp);247 __CLC_LONGN mantissa1 = temp_ux & 0x000FFFFFFFFFFFFFL;248 249 xexp = __CLC_CONVERT_LONGN(mask_exp_1023) ? xexp1 : xexp;250 mantissa = __CLC_CONVERT_LONGN(mask_exp_1023) ? mantissa1 : mantissa;251 252 __CLC_LONGN rax = (mantissa & 0x000ff00000000000) +253 ((mantissa & 0x0000080000000000) << 1);254 __CLC_INTN index = __CLC_CONVERT_INTN(rax >> 44);255 256 __CLC_GENTYPE F = __CLC_AS_GENTYPE(rax | 0x3FE0000000000000L);257 __CLC_GENTYPE Y = __CLC_AS_GENTYPE(mantissa | 0x3FE0000000000000L);258 __CLC_GENTYPE f = F - Y;259 __CLC_GENTYPE log_h = __CLC_USE_TABLE(log_f_inv_tbl_head, index);260 __CLC_GENTYPE log_t = __CLC_USE_TABLE(log_f_inv_tbl_tail, index);261 __CLC_GENTYPE f_inv = (log_h + log_t) * f;262 __CLC_GENTYPE r1 =263 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(f_inv) & 0xfffffffff8000000L);264 __CLC_GENTYPE r2 = __clc_fma(-F, r1, f) * (log_h + log_t);265 __CLC_GENTYPE r = r1 + r2;266 267 __CLC_GENTYPE poly = __clc_fma(268 r,269 __clc_fma(r,270 __clc_fma(r, __clc_fma(r, 1.0 / 7.0, 1.0 / 6.0), 1.0 / 5.0),271 1.0 / 4.0),272 1.0 / 3.0);273 poly = poly * r * r * r;274 275 __CLC_GENTYPE hr1r1 = 0.5 * r1 * r1;276 __CLC_GENTYPE poly0h = r1 + hr1r1;277 __CLC_GENTYPE poly0t = r1 - poly0h + hr1r1;278 poly = __clc_fma(r1, r2, __clc_fma(0.5 * r2, r2, poly)) + r2 + poly0t;279 280 log_h = __CLC_USE_TABLE(powlog_tbl_head, index);281 log_t = __CLC_USE_TABLE(powlog_tbl_tail, index);282 283 __CLC_GENTYPE resT_t = __clc_fma(xexp, real_log2_tail, +log_t) - poly;284 __CLC_GENTYPE resT = resT_t - poly0h;285 __CLC_GENTYPE resH = __clc_fma(xexp, real_log2_lead, log_h);286 __CLC_GENTYPE resT_h = poly0h;287 288 __CLC_GENTYPE H = resT + resH;289 __CLC_GENTYPE H_h =290 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(H) & 0xfffffffff8000000L);291 __CLC_GENTYPE T =292 (resH - H + resT) + (resT_t - (resT + resT_h)) + (H - H_h);293 H = H_h;294 295 __CLC_GENTYPE y_head =296 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(uy) & 0xfffffffff8000000L);297 __CLC_GENTYPE y_tail = y - y_head;298 299 __CLC_BIT_INTN mask_2_24 = ay > 0x4170000000000000; // 2^24300 __CLC_INTN nyh = __CLC_CONVERT_INTN(y_head);301 __CLC_INTN nyt = ny - nyh;302 __CLC_GENTYPE y_tail1 = __CLC_CONVERT_GENTYPE(nyt);303 y_tail = mask_2_24 ? y_tail1 : y_tail;304 305 __CLC_GENTYPE temp = __clc_fma(y_tail, H, __clc_fma(y_head, T, y_tail * T));306 v = __clc_fma(y_head, H, temp);307 vt = __clc_fma(y_head, H, -v) + temp;308 }309 310 // Now calculate exp of (v,vt)311 312 __CLC_GENTYPE expv;313 {314 const __CLC_GENTYPE max_exp_arg = 709.782712893384;315 const __CLC_GENTYPE min_exp_arg = -745.1332191019411;316 const __CLC_GENTYPE sixtyfour_by_lnof2 = 92.33248261689366;317 const __CLC_GENTYPE lnof2_by_64_head = 0.010830424260348081;318 const __CLC_GENTYPE lnof2_by_64_tail = -4.359010638708991e-10;319 320 // If v is so large that we need to return INFINITY, or so small that we321 // need to return 0, set v to known values that will produce that result. Do322 // not try to continue the computation with the original v and patch it up323 // afterwards because v may be so large that temp is out of range of int, in324 // which case that conversion, and a value based on that conversion being325 // passed to __clc_ldexp, results in undefined behavior.326 v = v > max_exp_arg ? 1000.0 : v;327 v = v < min_exp_arg ? -1000.0 : v;328 329 __CLC_GENTYPE temp = v * sixtyfour_by_lnof2;330 __CLC_INTN n = __CLC_CONVERT_INTN(temp);331 __CLC_GENTYPE dn = __CLC_CONVERT_GENTYPE(n);332 __CLC_INTN j = n & 0x0000003f;333 __CLC_INTN m = n >> 6;334 335 __CLC_GENTYPE f1 = __CLC_USE_TABLE(two_to_jby64_ep_tbl_head, j);336 __CLC_GENTYPE f2 = __CLC_USE_TABLE(two_to_jby64_ep_tbl_tail, j);337 __CLC_GENTYPE f = f1 + f2;338 339 __CLC_GENTYPE r1 = __clc_fma(dn, -lnof2_by_64_head, v);340 __CLC_GENTYPE r2 = dn * lnof2_by_64_tail;341 __CLC_GENTYPE r = (r1 + r2) + vt;342 343 __CLC_GENTYPE q =344 __clc_fma(r,345 __clc_fma(r,346 __clc_fma(r,347 __clc_fma(r, 1.38889490863777199667e-03,348 8.33336798434219616221e-03),349 4.16666666662260795726e-02),350 1.66666666665260878863e-01),351 5.00000000000000008883e-01);352 q = __clc_fma(r * r, q, r);353 354 expv = __clc_fma(f, q, f2) + f1;355 expv = __clc_ldexp(expv, m);356 }357 358 // See whether y is an integer.359 // inty = 0 means not an integer.360 // inty = 1 means odd integer.361 // inty = 2 means even integer.362 363 __CLC_LONGN inty = __CLC_CONVERT_LONGN(2 - (ny & 1));364 365 expv *= ((inty == 1) && !xpos) ? -1.0 : 1.0;366 367 __CLC_LONGN ret = __CLC_AS_LONGN(expv);368 369 // Now all the edge cases370 __CLC_BIT_INTN x_is_ninf = ux == (__CLC_LONGN)NINFBITPATT_DP64;371 __CLC_BIT_INTN x_is_pinf = ux == (__CLC_LONGN)PINFBITPATT_DP64;372 __CLC_LONGN xinf =373 xpos ? (__CLC_LONGN)PINFBITPATT_DP64 : (__CLC_LONGN)NINFBITPATT_DP64;374 375 ret = ((ax == 0L) && !ypos && (inty == 1)) ? xinf : ret;376 ret = ((ax == 0L) && !ypos && (inty == 2)) ? (__CLC_LONGN)PINFBITPATT_DP64377 : ret;378 ret = ((ax == 0L) && ypos && (inty == 2)) ? 0L : ret;379 __CLC_LONGN xzero = !xpos ? (__CLC_LONGN)0x8000000000000000L : 0L;380 ret = ((ax == 0L) && ypos && (inty == 1)) ? xzero : ret;381 ret = (x_is_ninf && !ypos && (inty == 1)) ? (__CLC_LONGN)0x8000000000000000L382 : ret;383 ret = (x_is_ninf && !ypos && (inty != 1)) ? 0L : ret;384 ret =385 (x_is_ninf && ypos && (inty == 1)) ? (__CLC_LONGN)NINFBITPATT_DP64 : ret;386 ret =387 (x_is_ninf && ypos && (inty != 1)) ? (__CLC_LONGN)PINFBITPATT_DP64 : ret;388 ret = (x_is_pinf && !ypos) ? 0L : ret;389 ret = (x_is_pinf && ypos) ? (__CLC_LONGN)PINFBITPATT_DP64 : ret;390 ret = ax > (__CLC_LONGN)PINFBITPATT_DP64 ? ux : ret;391 ret = __CLC_CONVERT_LONGN(ny == 0) ? (__CLC_LONGN)0x3ff0000000000000L : ret;392 393 return __CLC_AS_GENTYPE(ret);394}395 396#elif __CLC_FPSIZE == 16397 398_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_pown(__CLC_GENTYPE x, __CLC_INTN y) {399 return __CLC_CONVERT_GENTYPE(__clc_pown(__CLC_CONVERT_FLOATN(x), y));400}401 402#endif403