170 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/* Refer to the exp routine for the underlying algorithm */10#if __CLC_FPSIZE == 3211 12_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_expm1(__CLC_GENTYPE x) {13 // 128*log2 : 88.72283911167314 const __CLC_GENTYPE X_MAX = 0x1.62e42ep+6f;15 // -149*log2 : -103.2789299034318416 const __CLC_GENTYPE X_MIN = -0x1.9d1da0p+6f;17 // 64/log2 : 92.33248261689365718 const __CLC_GENTYPE R_64_BY_LOG2 = 0x1.715476p+6f;19 // log2/64 lead: 0.010803222720 const __CLC_GENTYPE R_LOG2_BY_64_LD = 0x1.620000p-7f;21 // log2/64 tail: 0.000027202038822 const __CLC_GENTYPE R_LOG2_BY_64_TL = 0x1.c85fdep-16f;23 24 __CLC_UINTN xi = __CLC_AS_UINTN(x);25 __CLC_INTN n = __CLC_CONVERT_INTN(x * R_64_BY_LOG2);26 __CLC_GENTYPE fn = __CLC_CONVERT_GENTYPE(n);27 28 __CLC_INTN j = n & 0x3f;29 __CLC_INTN m = n >> 6;30 31 __CLC_GENTYPE r =32 __clc_mad(fn, -R_LOG2_BY_64_TL, __clc_mad(fn, -R_LOG2_BY_64_LD, x));33 34 // Truncated Taylor series35 __CLC_GENTYPE z2 = __clc_mad(36 r * r, __clc_mad(r, __clc_mad(r, 0x1.555556p-5f, 0x1.555556p-3f), 0.5f),37 r);38 39 __CLC_GENTYPE m2 = __CLC_AS_GENTYPE((m + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);40 __CLC_GENTYPE exp_head = __CLC_USE_TABLE(exp_tbl_ep_head, j);41 __CLC_GENTYPE exp_tail = __CLC_USE_TABLE(exp_tbl_ep_tail, j);42 43 __CLC_GENTYPE two_to_jby64_h = exp_head * m2;44 __CLC_GENTYPE two_to_jby64_t = exp_tail * m2;45 __CLC_GENTYPE two_to_jby64 = two_to_jby64_h + two_to_jby64_t;46 47 z2 = __clc_mad(z2, two_to_jby64, two_to_jby64_t) + (two_to_jby64_h - 1.0f);48 // Make subnormals work49 z2 = x == 0.f ? x : z2;50 z2 = x < X_MIN || m < -24 ? -1.0f : z2;51 z2 = x > X_MAX ? __CLC_AS_GENTYPE((__CLC_UINTN)PINFBITPATT_SP32) : z2;52 z2 = __clc_isnan(x) ? x : z2;53 54 return z2;55}56 57#elif __CLC_FPSIZE == 6458 59_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_expm1(__CLC_GENTYPE x) {60 const __CLC_GENTYPE max_expm1_arg = 709.8;61 const __CLC_GENTYPE min_expm1_arg = -37.42994775023704;62 // 0x3FCC8FF7C79A9A22 = log(1+1/4)63 const __CLC_GENTYPE log_OnePlus_OneByFour = 0.22314355131420976;64 // 0xBFD269621134DB93 = log(1-1/4)65 const __CLC_GENTYPE log_OneMinus_OneByFour = -0.28768207245178096;66 const __CLC_GENTYPE sixtyfour_by_lnof2 =67 92.33248261689366; // 0x40571547652b82fe68 const __CLC_GENTYPE lnof2_by_64_head =69 0.010830424696223417; // 0x3f862e42fefa000070 const __CLC_GENTYPE lnof2_by_64_tail =71 2.5728046223276688e-14; // 0x3d1cf79abc9e3b3972 73 // First, assume log(1-1/4) < x < log(1+1/4) i.e -0.28768 < x < 0.2231474 __CLC_GENTYPE u = __CLC_AS_GENTYPE(__CLC_AS_ULONGN(x) & 0xffffffffff000000UL);75 __CLC_GENTYPE v = x - u;76 __CLC_GENTYPE y = u * u * 0.5;77 __CLC_GENTYPE z = v * (x + u) * 0.5;78 79 __CLC_GENTYPE q = __clc_fma(80 x,81 __clc_fma(82 x,83 __clc_fma(84 x,85 __clc_fma(86 x,87 __clc_fma(88 x,89 __clc_fma(x,90 __clc_fma(x,91 __clc_fma(x, 2.4360682937111612e-8,92 2.7582184028154370e-7),93 2.7558212415361945e-6),94 2.4801576918453420e-5),95 1.9841269447671544e-4),96 1.3888888890687830e-3),97 8.3333333334012270e-3),98 4.1666666666665560e-2),99 1.6666666666666632e-1);100 q *= x * x * x;101 102 __CLC_GENTYPE z1g = (u + y) + (q + (v + z));103 __CLC_GENTYPE z1 = x + (y + (q + z));104 z1 = y >= 0x1.0p-7 ? z1g : z1;105 106 // Now assume outside interval around 0107 __CLC_INTN n = __CLC_CONVERT_INTN(x * sixtyfour_by_lnof2);108 __CLC_INTN j = n & 0x3f;109 __CLC_INTN m = n >> 6;110 111 __CLC_GENTYPE f1 = __CLC_USE_TABLE(two_to_jby64_ep_tbl_head, j);112 __CLC_GENTYPE f2 = __CLC_USE_TABLE(two_to_jby64_ep_tbl_tail, j);113 __CLC_GENTYPE f = f1 + f2;114 115 __CLC_GENTYPE dn = __CLC_CONVERT_GENTYPE(-n);116 __CLC_GENTYPE r =117 __clc_fma(dn, lnof2_by_64_tail, __clc_fma(dn, lnof2_by_64_head, x));118 119 q = __clc_fma(r,120 __clc_fma(r,121 __clc_fma(r,122 __clc_fma(r, 1.38889490863777199667e-03,123 8.33336798434219616221e-03),124 4.16666666662260795726e-02),125 1.66666666665260878863e-01),126 5.00000000000000008883e-01);127 q = __clc_fma(r * r, q, r);128 129 __CLC_GENTYPE twopm = __CLC_AS_GENTYPE(__CLC_CONVERT_LONGN(m + EXPBIAS_DP64)130 << EXPSHIFTBITS_DP64);131 __CLC_GENTYPE twopmm = __CLC_AS_GENTYPE(__CLC_CONVERT_LONGN(EXPBIAS_DP64 - m)132 << EXPSHIFTBITS_DP64);133 134 // Computations for m > 52, including where result is close to Inf135 __CLC_ULONGN uval = __CLC_AS_ULONGN(0x1.0p+1023 * (f1 + (f * q + (f2))));136 __CLC_INTN e = __CLC_CONVERT_INTN(uval >> EXPSHIFTBITS_DP64) + 1;137 138 __CLC_GENTYPE zme1024 = __CLC_AS_GENTYPE(139 (__CLC_CONVERT_ULONGN(e) << EXPSHIFTBITS_DP64) | (uval & MANTBITS_DP64));140 zme1024 = __CLC_CONVERT_LONGN(e == 2047)141 ? __CLC_AS_GENTYPE((__CLC_ULONGN)PINFBITPATT_DP64)142 : zme1024;143 144 __CLC_GENTYPE zmg52 = twopm * (f1 + __clc_fma(f, q, f2 - twopmm));145 zmg52 = __CLC_CONVERT_LONGN(m == 1024) ? zme1024 : zmg52;146 147 // For m < 53148 __CLC_GENTYPE zml53 =149 twopm * ((f1 - twopmm) + __clc_fma(f1, q, f2 * (1.0 + q)));150 151 // For m < -7152 __CLC_GENTYPE zmln7 = __clc_fma(twopm, f1 + __clc_fma(f, q, f2), -1.0);153 154 z = __CLC_CONVERT_LONGN(m < 53) ? zml53 : zmg52;155 z = __CLC_CONVERT_LONGN(m < -7) ? zmln7 : z;156 z = x > log_OneMinus_OneByFour && x < log_OnePlus_OneByFour ? z1 : z;157 z = x > max_expm1_arg ? __CLC_AS_GENTYPE((__CLC_ULONGN)PINFBITPATT_DP64) : z;158 z = x < min_expm1_arg ? -1.0 : z;159 160 return z;161}162 163#elif __CLC_FPSIZE == 16164 165_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_expm1(__CLC_GENTYPE x) {166 return __CLC_CONVERT_GENTYPE(__clc_expm1(__CLC_CONVERT_FLOATN(x)));167}168 169#endif170