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/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 // This is an algorithm for log10(x) in single precision which is 19 // correctly rounded for all rounding modes. 20 // - An exhaustive test show that when x >= 2^45, log1pf(x) == logf(x) 21 // for all rounding modes. 22 // - When 2^(-8) <= |x| < 2^45, the sum (double(x) + 1.0) is exact, 23 // so we can adapt the correctly rounded algorithm of logf to compute 24 // log(double(x) + 1.0) correctly. For more information about the logf 25 // algorithm, see `libc/src/math/generic/logf.cpp`. 26 // - When |x| < 2^(-8), we use a degree-6 polynomial in double precision 27 // generated with Sollya using the following command: 28 // fpminimax(log(1 + x)/x, 5, [|D...|], [-2^-8; 2^-8]); 29 30 namespace __llvm_libc { 31 32 namespace internal { 33 34 // We don't need to treat denormal 35 INLINE_FMA static inline float log(double x) { 36 constexpr double LOG_2 = 0x1.62e42fefa39efp-1; 37 38 using FPBits = typename fputil::FPBits<double>; 39 FPBits xbits(x); 40 41 if (xbits.is_zero()) { 42 return static_cast<float>(fputil::FPBits<float>::neg_inf()); 43 } 44 45 if (xbits.uintval() > FPBits::MAX_NORMAL) { 46 if (xbits.get_sign() && !xbits.is_nan()) { 47 return fputil::FPBits<float>::build_nan( 48 1 << (fputil::MantissaWidth<float>::VALUE - 1)); 49 } 50 return static_cast<float>(x); 51 } 52 53 double m = static_cast<double>(xbits.get_exponent()); 54 55 // Set bits to 1.m 56 xbits.set_unbiased_exponent(0x3FF); 57 // Get the 8 highest bits, use 7 bits (excluding the implicit hidden bit) for 58 // lookup tables. 59 int f_index = 60 xbits.get_mantissa() >> 45; // fputil::MantissaWidth<double>::VALUE - 7 61 62 FPBits f = xbits; 63 // Clear the lowest 45 bits. 64 f.bits &= ~0x0000'1FFF'FFFF'FFFFULL; 65 66 double d = static_cast<double>(xbits) - static_cast<double>(f); 67 d *= ONE_OVER_F[f_index]; 68 69 double extra_factor = fputil::multiply_add(m, LOG_2, LOG_F[f_index]); 70 71 double r = fputil::polyeval(d, extra_factor, 0x1.fffffffffffacp-1, 72 -0x1.fffffffef9cb2p-2, 0x1.5555513bc679ap-2, 73 -0x1.fff4805ea441p-3, 0x1.930180dbde91ap-3); 74 75 return static_cast<float>(r); 76 } 77 78 } // namespace internal 79 80 INLINE_FMA 81 LLVM_LIBC_FUNCTION(float, log1pf, (float x)) { 82 using FPBits = typename fputil::FPBits<float>; 83 FPBits xbits(x); 84 double xd = static_cast<double>(x); 85 86 if (xbits.get_exponent() >= -8) { 87 // Hard-to-round cases. 88 switch (xbits.uintval()) { 89 case 0x3b9315c8U: // x = 0x1.262b9p-8f 90 if (fputil::get_round() != FE_UPWARD) 91 return 0x1.25830cp-8f; 92 break; 93 case 0x3c6eb7afU: // x = 0x1.dd6f5ep-7f 94 if (fputil::get_round() == FE_UPWARD) 95 return 0x1.d9fd86p-7f; 96 return 0x1.d9fd84p-7f; 97 case 0x41078febU: // x = 0x1.0f1fd6p+3f 98 if (fputil::get_round() != FE_UPWARD) 99 return 0x1.1fcbcep+1f; 100 break; 101 case 0x5cd69e88U: // x = 0x1.ad3d1p+58f 102 if (fputil::get_round() != FE_UPWARD) 103 return 0x1.45c146p+5f; 104 break; 105 case 0x65d890d3U: // x = 0x1.b121a6p+76f 106 if (fputil::get_round() == FE_TONEAREST) 107 return 0x1.a9a3f2p+5f; 108 break; 109 case 0x6f31a8ecU: // x = 0x1.6351d8p+95f 110 if (fputil::get_round() == FE_TONEAREST) 111 return 0x1.08b512p+6f; 112 break; 113 case 0x7a17f30aU: // x = 0x1.2fe614p+117f 114 if (fputil::get_round() != FE_UPWARD) 115 return 0x1.451436p+6f; 116 break; 117 case 0xbc4d092cU: // x = -0x1.9a1258p-7f 118 if (fputil::get_round() == FE_TONEAREST) 119 return -0x1.9ca8bep-7f; 120 break; 121 case 0xbc657728U: // x = -0x1.caee5p-7f 122 if (fputil::get_round() != FE_DOWNWARD) 123 return -0x1.ce2cccp-7f; 124 break; 125 case 0xbd1d20afU: // x = -0x1.3a415ep-5f 126 int round_mode = fputil::get_round(); 127 if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO) 128 return -0x1.40711p-5f; 129 return -0x1.407112p-5f; 130 } 131 132 return internal::log(xd + 1.0); 133 } 134 135 // Hard-to round cases. 136 switch (xbits.uintval()) { 137 case 0x35400003U: // x = 0x1.800006p-21f 138 if (fputil::get_round() == FE_TONEAREST) 139 return 0x1.7ffffep-21f; 140 break; 141 case 0x3710001bU: // x = 0x1.200036p-17f 142 if (fputil::get_round() == FE_TONEAREST) 143 return 0x1.1fffe6p-17f; 144 break; 145 case 0xb53ffffdU: // x = -0x1.7ffffap-21f 146 if (fputil::get_round() != FE_DOWNWARD) 147 return -0x1.800002p-21f; 148 break; 149 case 0xb70fffe5U: // x = -0x1.1fffcap-17f 150 if (fputil::get_round() != FE_DOWNWARD) 151 return -0x1.20001ap-17f; 152 break; 153 case 0xbb0ec8c4U: // x = -0x1.1d9188p-9f 154 if (fputil::get_round() == FE_TONEAREST) 155 return -0x1.1de14ap-9f; 156 break; 157 } 158 159 double r; 160 // Polymial generated with Sollya: 161 // > fpminimax(log(1 + x)/x, 5, [|D...|], [-2^-8; 2^-8]); 162 r = fputil::polyeval(xd, -0x1p-1, 0x1.5555555515551p-2, -0x1.ffffffff82bdap-3, 163 0x1.999b33348d3aep-3, -0x1.5556cae3adcc3p-3); 164 return static_cast<float>(fputil::multiply_add(r, xd * xd, xd)); 165 } 166 167 } // namespace __llvm_libc 168