//===----------------------------------------------------------------------===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// /* Refer to the exp routine for the underlying algorithm */ #if __CLC_FPSIZE == 32 _CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_expm1(__CLC_GENTYPE x) { // 128*log2 : 88.722839111673 const __CLC_GENTYPE X_MAX = 0x1.62e42ep+6f; // -149*log2 : -103.27892990343184 const __CLC_GENTYPE X_MIN = -0x1.9d1da0p+6f; // 64/log2 : 92.332482616893657 const __CLC_GENTYPE R_64_BY_LOG2 = 0x1.715476p+6f; // log2/64 lead: 0.0108032227 const __CLC_GENTYPE R_LOG2_BY_64_LD = 0x1.620000p-7f; // log2/64 tail: 0.0000272020388 const __CLC_GENTYPE R_LOG2_BY_64_TL = 0x1.c85fdep-16f; __CLC_UINTN xi = __CLC_AS_UINTN(x); __CLC_INTN n = __CLC_CONVERT_INTN(x * R_64_BY_LOG2); __CLC_GENTYPE fn = __CLC_CONVERT_GENTYPE(n); __CLC_INTN j = n & 0x3f; __CLC_INTN m = n >> 6; __CLC_GENTYPE r = __clc_mad(fn, -R_LOG2_BY_64_TL, __clc_mad(fn, -R_LOG2_BY_64_LD, x)); // Truncated Taylor series __CLC_GENTYPE z2 = __clc_mad( r * r, __clc_mad(r, __clc_mad(r, 0x1.555556p-5f, 0x1.555556p-3f), 0.5f), r); __CLC_GENTYPE m2 = __CLC_AS_GENTYPE((m + EXPBIAS_SP32) << EXPSHIFTBITS_SP32); __CLC_GENTYPE exp_head = __CLC_USE_TABLE(exp_tbl_ep_head, j); __CLC_GENTYPE exp_tail = __CLC_USE_TABLE(exp_tbl_ep_tail, j); __CLC_GENTYPE two_to_jby64_h = exp_head * m2; __CLC_GENTYPE two_to_jby64_t = exp_tail * m2; __CLC_GENTYPE two_to_jby64 = two_to_jby64_h + two_to_jby64_t; z2 = __clc_mad(z2, two_to_jby64, two_to_jby64_t) + (two_to_jby64_h - 1.0f); // Make subnormals work z2 = x == 0.f ? x : z2; z2 = x < X_MIN || m < -24 ? -1.0f : z2; z2 = x > X_MAX ? __CLC_AS_GENTYPE((__CLC_UINTN)PINFBITPATT_SP32) : z2; z2 = __clc_isnan(x) ? x : z2; return z2; } #elif __CLC_FPSIZE == 64 _CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_expm1(__CLC_GENTYPE x) { const __CLC_GENTYPE max_expm1_arg = 709.8; const __CLC_GENTYPE min_expm1_arg = -37.42994775023704; // 0x3FCC8FF7C79A9A22 = log(1+1/4) const __CLC_GENTYPE log_OnePlus_OneByFour = 0.22314355131420976; // 0xBFD269621134DB93 = log(1-1/4) const __CLC_GENTYPE log_OneMinus_OneByFour = -0.28768207245178096; const __CLC_GENTYPE sixtyfour_by_lnof2 = 92.33248261689366; // 0x40571547652b82fe const __CLC_GENTYPE lnof2_by_64_head = 0.010830424696223417; // 0x3f862e42fefa0000 const __CLC_GENTYPE lnof2_by_64_tail = 2.5728046223276688e-14; // 0x3d1cf79abc9e3b39 // First, assume log(1-1/4) < x < log(1+1/4) i.e -0.28768 < x < 0.22314 __CLC_GENTYPE u = __CLC_AS_GENTYPE(__CLC_AS_ULONGN(x) & 0xffffffffff000000UL); __CLC_GENTYPE v = x - u; __CLC_GENTYPE y = u * u * 0.5; __CLC_GENTYPE z = v * (x + u) * 0.5; __CLC_GENTYPE q = __clc_fma( x, __clc_fma( x, __clc_fma( x, __clc_fma( x, __clc_fma( x, __clc_fma(x, __clc_fma(x, __clc_fma(x, 2.4360682937111612e-8, 2.7582184028154370e-7), 2.7558212415361945e-6), 2.4801576918453420e-5), 1.9841269447671544e-4), 1.3888888890687830e-3), 8.3333333334012270e-3), 4.1666666666665560e-2), 1.6666666666666632e-1); q *= x * x * x; __CLC_GENTYPE z1g = (u + y) + (q + (v + z)); __CLC_GENTYPE z1 = x + (y + (q + z)); z1 = y >= 0x1.0p-7 ? z1g : z1; // Now assume outside interval around 0 __CLC_INTN n = __CLC_CONVERT_INTN(x * sixtyfour_by_lnof2); __CLC_INTN j = n & 0x3f; __CLC_INTN m = n >> 6; __CLC_GENTYPE f1 = __CLC_USE_TABLE(two_to_jby64_ep_tbl_head, j); __CLC_GENTYPE f2 = __CLC_USE_TABLE(two_to_jby64_ep_tbl_tail, j); __CLC_GENTYPE f = f1 + f2; __CLC_GENTYPE dn = __CLC_CONVERT_GENTYPE(-n); __CLC_GENTYPE r = __clc_fma(dn, lnof2_by_64_tail, __clc_fma(dn, lnof2_by_64_head, x)); q = __clc_fma(r, __clc_fma(r, __clc_fma(r, __clc_fma(r, 1.38889490863777199667e-03, 8.33336798434219616221e-03), 4.16666666662260795726e-02), 1.66666666665260878863e-01), 5.00000000000000008883e-01); q = __clc_fma(r * r, q, r); __CLC_GENTYPE twopm = __CLC_AS_GENTYPE(__CLC_CONVERT_LONGN(m + EXPBIAS_DP64) << EXPSHIFTBITS_DP64); __CLC_GENTYPE twopmm = __CLC_AS_GENTYPE(__CLC_CONVERT_LONGN(EXPBIAS_DP64 - m) << EXPSHIFTBITS_DP64); // Computations for m > 52, including where result is close to Inf __CLC_ULONGN uval = __CLC_AS_ULONGN(0x1.0p+1023 * (f1 + (f * q + (f2)))); __CLC_INTN e = __CLC_CONVERT_INTN(uval >> EXPSHIFTBITS_DP64) + 1; __CLC_GENTYPE zme1024 = __CLC_AS_GENTYPE( (__CLC_CONVERT_ULONGN(e) << EXPSHIFTBITS_DP64) | (uval & MANTBITS_DP64)); zme1024 = __CLC_CONVERT_LONGN(e == 2047) ? __CLC_AS_GENTYPE((__CLC_ULONGN)PINFBITPATT_DP64) : zme1024; __CLC_GENTYPE zmg52 = twopm * (f1 + __clc_fma(f, q, f2 - twopmm)); zmg52 = __CLC_CONVERT_LONGN(m == 1024) ? zme1024 : zmg52; // For m < 53 __CLC_GENTYPE zml53 = twopm * ((f1 - twopmm) + __clc_fma(f1, q, f2 * (1.0 + q))); // For m < -7 __CLC_GENTYPE zmln7 = __clc_fma(twopm, f1 + __clc_fma(f, q, f2), -1.0); z = __CLC_CONVERT_LONGN(m < 53) ? zml53 : zmg52; z = __CLC_CONVERT_LONGN(m < -7) ? zmln7 : z; z = x > log_OneMinus_OneByFour && x < log_OnePlus_OneByFour ? z1 : z; z = x > max_expm1_arg ? __CLC_AS_GENTYPE((__CLC_ULONGN)PINFBITPATT_DP64) : z; z = x < min_expm1_arg ? -1.0 : z; return z; } #elif __CLC_FPSIZE == 16 _CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_expm1(__CLC_GENTYPE x) { return __CLC_CONVERT_GENTYPE(__clc_expm1(__CLC_CONVERT_FLOATN(x))); } #endif