1 //===-- Single-precision log1p(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/log1pf.h" 10 #include "common_constants.h" // Lookup table for (1/f) and log(f) 11 #include "src/__support/FPUtil/FEnvImpl.h" 12 #include "src/__support/FPUtil/FMA.h" 13 #include "src/__support/FPUtil/FPBits.h" 14 #include "src/__support/FPUtil/PolyEval.h" 15 #include "src/__support/FPUtil/except_value_utils.h" 16 #include "src/__support/FPUtil/multiply_add.h" 17 #include "src/__support/common.h" 18 #include "src/__support/macros/config.h" 19 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY 20 #include "src/__support/macros/properties/cpu_features.h" 21 22 // This is an algorithm for log10(x) in single precision which is 23 // correctly rounded for all rounding modes. 24 // - An exhaustive test show that when x >= 2^45, log1pf(x) == logf(x) 25 // for all rounding modes. 26 // - When 2^(-6) <= |x| < 2^45, the sum (double(x) + 1.0) is exact, 27 // so we can adapt the correctly rounded algorithm of logf to compute 28 // log(double(x) + 1.0) correctly. For more information about the logf 29 // algorithm, see `libc/src/math/generic/logf.cpp`. 30 // - When |x| < 2^(-6), we use a degree-8 polynomial in double precision 31 // generated with Sollya using the following command: 32 // fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); 33 34 namespace LIBC_NAMESPACE_DECL { 35 36 namespace internal { 37 38 // We don't need to treat denormal and 0 39 LIBC_INLINE float log(double x) { 40 constexpr double LOG_2 = 0x1.62e42fefa39efp-1; 41 42 using FPBits = typename fputil::FPBits<double>; 43 FPBits xbits(x); 44 45 uint64_t x_u = xbits.uintval(); 46 47 if (LIBC_UNLIKELY(x_u > FPBits::max_normal().uintval())) { 48 if (xbits.is_neg() && !xbits.is_nan()) { 49 fputil::set_errno_if_required(EDOM); 50 fputil::raise_except_if_required(FE_INVALID); 51 return fputil::FPBits<float>::quiet_nan().get_val(); 52 } 53 return static_cast<float>(x); 54 } 55 56 double m = static_cast<double>(xbits.get_exponent()); 57 58 // Get the 8 highest bits, use 7 bits (excluding the implicit hidden bit) for 59 // lookup tables. 60 int f_index = static_cast<int>(xbits.get_mantissa() >> 61 (fputil::FPBits<double>::FRACTION_LEN - 7)); 62 63 // Set bits to 1.m 64 xbits.set_biased_exponent(0x3FF); 65 FPBits f = xbits; 66 67 // Clear the lowest 45 bits. 68 f.set_uintval(f.uintval() & ~0x0000'1FFF'FFFF'FFFFULL); 69 70 double d = xbits.get_val() - f.get_val(); 71 d *= ONE_OVER_F[f_index]; 72 73 double extra_factor = fputil::multiply_add(m, LOG_2, LOG_F[f_index]); 74 75 double r = fputil::polyeval(d, extra_factor, 0x1.fffffffffffacp-1, 76 -0x1.fffffffef9cb2p-2, 0x1.5555513bc679ap-2, 77 -0x1.fff4805ea441p-3, 0x1.930180dbde91ap-3); 78 79 return static_cast<float>(r); 80 } 81 82 } // namespace internal 83 84 LLVM_LIBC_FUNCTION(float, log1pf, (float x)) { 85 using FPBits = typename fputil::FPBits<float>; 86 FPBits xbits(x); 87 uint32_t x_u = xbits.uintval(); 88 uint32_t x_a = x_u & 0x7fff'ffffU; 89 double xd = static_cast<double>(x); 90 91 // Use log1p(x) = log(1 + x) for |x| > 2^-6; 92 if (x_a > 0x3c80'0000U) { 93 // Hard-to-round cases. 94 switch (x_u) { 95 case 0x41078febU: // x = 0x1.0f1fd6p3 96 return fputil::round_result_slightly_up(0x1.1fcbcep1f); 97 case 0x5cd69e88U: // x = 0x1.ad3d1p+58f 98 return fputil::round_result_slightly_up(0x1.45c146p+5f); 99 case 0x65d890d3U: // x = 0x1.b121a6p+76f 100 return fputil::round_result_slightly_down(0x1.a9a3f2p+5f); 101 case 0x6f31a8ecU: // x = 0x1.6351d8p+95f 102 return fputil::round_result_slightly_down(0x1.08b512p+6f); 103 case 0x7a17f30aU: // x = 0x1.2fe614p+117f 104 return fputil::round_result_slightly_up(0x1.451436p+6f); 105 case 0xbd1d20afU: // x = -0x1.3a415ep-5f 106 return fputil::round_result_slightly_up(-0x1.407112p-5f); 107 case 0xbf800000U: // x = -1.0 108 fputil::set_errno_if_required(ERANGE); 109 fputil::raise_except_if_required(FE_DIVBYZERO); 110 return FPBits::inf(Sign::NEG).get_val(); 111 #ifndef LIBC_TARGET_CPU_HAS_FMA 112 case 0x4cc1c80bU: // x = 0x1.839016p+26f 113 return fputil::round_result_slightly_down(0x1.26fc04p+4f); 114 case 0x5ee8984eU: // x = 0x1.d1309cp+62f 115 return fputil::round_result_slightly_up(0x1.5c9442p+5f); 116 case 0x665e7ca6U: // x = 0x1.bcf94cp+77f 117 return fputil::round_result_slightly_up(0x1.af66cp+5f); 118 case 0x79e7ec37U: // x = 0x1.cfd86ep+116f 119 return fputil::round_result_slightly_up(0x1.43ff6ep+6f); 120 #endif // LIBC_TARGET_CPU_HAS_FMA 121 } 122 123 return internal::log(xd + 1.0); 124 } 125 126 // |x| <= 2^-6. 127 // Hard-to round cases. 128 switch (x_u) { 129 case 0x35400003U: // x = 0x1.800006p-21f 130 return fputil::round_result_slightly_down(0x1.7ffffep-21f); 131 case 0x3710001bU: // x = 0x1.200036p-17f 132 return fputil::round_result_slightly_down(0x1.1fffe6p-17f); 133 case 0xb53ffffdU: // x = -0x1.7ffffap-21 134 return fputil::round_result_slightly_down(-0x1.800002p-21f); 135 case 0xb70fffe5U: // x = -0x1.1fffcap-17 136 return fputil::round_result_slightly_down(-0x1.20001ap-17f); 137 case 0xbb0ec8c4U: // x = -0x1.1d9188p-9 138 return fputil::round_result_slightly_up(-0x1.1de14ap-9f); 139 } 140 141 // Polymial generated by Sollya with: 142 // > fpminimax(log(1 + x)/x, 7, [|D...|], [-2^-6; 2^-6]); 143 const double COEFFS[7] = {-0x1.0000000000000p-1, 0x1.5555555556aadp-2, 144 -0x1.000000000181ap-2, 0x1.999998998124ep-3, 145 -0x1.55555452e2a2bp-3, 0x1.24adb8cde4aa7p-3, 146 -0x1.0019db915ef6fp-3}; 147 148 double xsq = xd * xd; 149 double c0 = fputil::multiply_add(xd, COEFFS[1], COEFFS[0]); 150 double c1 = fputil::multiply_add(xd, COEFFS[3], COEFFS[2]); 151 double c2 = fputil::multiply_add(xd, COEFFS[5], COEFFS[4]); 152 double r = fputil::polyeval(xsq, xd, c0, c1, c2, COEFFS[6]); 153 154 return static_cast<float>(r); 155 } 156 157 } // namespace LIBC_NAMESPACE_DECL 158