1 //===-- Single-precision e^x 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/expf.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, expf, (float x)) { 24 using FPBits = typename fputil::FPBits<float>; 25 FPBits xbits(x); 26 27 // When x < log(2^-150) or nan 28 if (unlikely(xbits.uintval() >= 0xc2cf'f1b5U)) { 29 // exp(-Inf) = 0 30 if (xbits.is_inf()) 31 return 0.0f; 32 // exp(nan) = nan 33 if (xbits.is_nan()) 34 return x; 35 if (fputil::get_round() == FE_UPWARD) 36 return static_cast<float>(FPBits(FPBits::MIN_SUBNORMAL)); 37 errno = ERANGE; 38 return 0.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 // |x| < 2^-25 52 if (unlikely(xbits.get_unbiased_exponent() <= 101)) { 53 return 1.0f + x; 54 } 55 56 // For -104 < x < 89, to compute exp(x), we perform the following range 57 // reduction: find hi, mid, lo such that: 58 // x = hi + mid + lo, in which 59 // hi is an integer, 60 // mid * 2^7 is an integer 61 // -2^(-8) <= lo < 2^-8. 62 // In particular, 63 // hi + mid = round(x * 2^7) * 2^(-7). 64 // Then, 65 // exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo). 66 // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2 67 // respectively. exp(lo) is computed using a degree-7 minimax polynomial 68 // generated by Sollya. 69 70 // x_hi = hi + mid. 71 int x_hi = static_cast<int>(x * 0x1.0p7f); 72 // Subtract (hi + mid) from x to get lo. 73 x -= static_cast<float>(x_hi) * 0x1.0p-7f; 74 double xd = static_cast<double>(x); 75 // Make sure that -2^(-8) <= lo < 2^-8. 76 if (x >= 0x1.0p-8f) { 77 ++x_hi; 78 xd -= 0x1.0p-7; 79 } 80 if (x < -0x1.0p-8f) { 81 --x_hi; 82 xd += 0x1.0p-7; 83 } 84 x_hi += 104 << 7; 85 // hi = x_hi >> 7 86 double exp_hi = EXP_M1[x_hi >> 7]; 87 // lo = x_hi & 0x0000'007fU; 88 double exp_mid = EXP_M2[x_hi & 0x7f]; 89 // Degree-7 minimax polynomial generated by Sollya with the following 90 // commands: 91 // > display = hexadecimal; 92 // > Q = fpminimax(expm1(x)/x, 6, [|D...|], [-2^-8, 2^-8]); 93 // > Q; 94 double exp_lo = fputil::polyeval( 95 xd, 0x1p0, 0x1p0, 0x1p-1, 0x1.5555555555555p-3, 0x1.55555555553ap-5, 96 0x1.1111111204dfcp-7, 0x1.6c16cb2da593ap-10, 0x1.9ff1648996d2ep-13); 97 return static_cast<float>(exp_hi * exp_mid * exp_lo); 98 } 99 100 } // namespace __llvm_libc 101