1 //===-- Single-precision e^x - 1 function ---------------------------------===// 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-exception 6 // 7 //===----------------------------------------------------------------------===// 8 9 #include "src/math/expm1f.h" 10 #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2. 11 #include "src/__support/FPUtil/BasicOperations.h" 12 #include "src/__support/FPUtil/FEnvImpl.h" 13 #include "src/__support/FPUtil/FMA.h" 14 #include "src/__support/FPUtil/FPBits.h" 15 #include "src/__support/FPUtil/PolyEval.h" 16 #include "src/__support/common.h" 17 18 #include <errno.h> 19 20 namespace __llvm_libc { 21 22 INLINE_FMA 23 LLVM_LIBC_FUNCTION(float, expm1f, (float x)) { 24 using FPBits = typename fputil::FPBits<float>; 25 FPBits xbits(x); 26 27 // When x < log(2^-25) or nan 28 if (unlikely(xbits.uintval() >= 0xc18a'a123U)) { 29 // exp(-Inf) = 0 30 if (xbits.is_inf()) 31 return -1.0f; 32 // exp(nan) = nan 33 if (xbits.is_nan()) 34 return x; 35 int round_mode = fputil::get_round(); 36 if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO) 37 return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f 38 return -1.0f; 39 } 40 // x >= 89 or nan 41 if (unlikely(!xbits.get_sign() && (xbits.uintval() >= 0x42b2'0000))) { 42 if (xbits.uintval() < 0x7f80'0000U) { 43 int rounding = fputil::get_round(); 44 if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) 45 return static_cast<float>(FPBits(FPBits::MAX_NORMAL)); 46 47 errno = ERANGE; 48 } 49 return x + static_cast<float>(FPBits::inf()); 50 } 51 52 int unbiased_exponent = static_cast<int>(xbits.get_unbiased_exponent()); 53 // |x| < 2^-4 54 if (unbiased_exponent < 123) { 55 // |x| < 2^-25 56 if (unbiased_exponent < 102) { 57 // x = -0.0f 58 if (unlikely(xbits.uintval() == 0x8000'0000U)) 59 return x; 60 // When |x| < 2^-25, the relative error: 61 // |(e^x - 1) - x| / |x| < |x^2| / |x| = |x| < 2^-25 < epsilon(1)/2. 62 // So the correctly rounded values of expm1(x) are: 63 // = x + eps(x) if rounding mode = FE_UPWARD, 64 // or (rounding mode = FE_TOWARDZERO and x is negative), 65 // = x otherwise. 66 // To simplify the rounding decision and make it more efficient, we use 67 // fma(x, x, x) ~ x + x^2 instead. 68 return fputil::fma(x, x, x); 69 } 70 // 2^-25 <= |x| < 2^-4 71 double xd = static_cast<double>(x); 72 double xsq = xd * xd; 73 // Degree-8 minimax polynomial generated by Sollya with: 74 // > display = hexadecimal; 75 // > P = fpminimax(expm1(x)/x, 7, [|D...|], [-2^-4, 2^-4]); 76 double r = 77 fputil::polyeval(xd, 0x1p-1, 0x1.55555555559abp-3, 0x1.55555555551a7p-5, 78 0x1.111110f70f2a4p-7, 0x1.6c16c17639e82p-10, 79 0x1.a02526febbea6p-13, 0x1.a01dc40888fcdp-16); 80 return static_cast<float>(fputil::fma(r, xsq, xd)); 81 } 82 83 // For -18 < x < 89, to compute exp(x), we perform the following range 84 // reduction: find hi, mid, lo such that: 85 // x = hi + mid + lo, in which 86 // hi is an integer, 87 // mid * 2^7 is an integer 88 // -2^(-8) <= lo < 2^-8. 89 // In particular, 90 // hi + mid = round(x * 2^7) * 2^(-7). 91 // Then, 92 // exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo). 93 // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2 94 // respectively. exp(lo) is computed using a degree-7 minimax polynomial 95 // generated by Sollya. 96 97 // Exceptional value 98 if (xbits.uintval() == 0xbdc1'c6cbU) { 99 // x = -0x1.838d96p-4f 100 int round_mode = fputil::get_round(); 101 if (round_mode == FE_TONEAREST || round_mode == FE_DOWNWARD) 102 return -0x1.71c884p-4f; 103 return -0x1.71c882p-4f; 104 } 105 106 // x_hi = hi + mid. 107 int x_hi = static_cast<int>(x * 0x1.0p7f); 108 // Subtract (hi + mid) from x to get lo. 109 x -= static_cast<float>(x_hi) * 0x1.0p-7f; 110 double xd = static_cast<double>(x); 111 // Make sure that -2^(-8) <= lo < 2^-8. 112 if (x >= 0x1.0p-8f) { 113 ++x_hi; 114 xd -= 0x1.0p-7; 115 } 116 if (x < -0x1.0p-8f) { 117 --x_hi; 118 xd += 0x1.0p-7; 119 } 120 x_hi += 104 << 7; 121 // hi = x_hi >> 7 122 double exp_hi = EXP_M1[x_hi >> 7]; 123 // lo = x_hi & 0x0000'007fU; 124 double exp_mid = EXP_M2[x_hi & 0x7f]; 125 double exp_hi_mid = exp_hi * exp_mid; 126 // Degree-7 minimax polynomial generated by Sollya with the following 127 // commands: 128 // > display = hexadecimal; 129 // > Q = fpminimax(expm1(x)/x, 6, [|D...|], [-2^-8, 2^-8]); 130 // > Q; 131 double exp_lo = fputil::polyeval( 132 xd, 0x1p0, 0x1p0, 0x1p-1, 0x1.5555555555555p-3, 0x1.55555555553ap-5, 133 0x1.1111111204dfcp-7, 0x1.6c16cb2da593ap-10, 0x1.9ff1648996d2ep-13); 134 return static_cast<float>(fputil::fma(exp_hi_mid, exp_lo, -1.0)); 135 } 136 137 } // namespace __llvm_libc 138