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 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 LLVM_LIBC_FUNCTION(float, log1pf, (float x)) { 81 using FPBits = typename fputil::FPBits<float>; 82 FPBits xbits(x); 83 double xd = static_cast<double>(x); 84 85 if (xbits.get_exponent() >= -8) { 86 // Hard-to-round cases. 87 switch (xbits.uintval()) { 88 case 0x3b9315c8U: // x = 0x1.262b9p-8f 89 if (fputil::get_round() != FE_UPWARD) 90 return 0x1.25830cp-8f; 91 break; 92 case 0x3c6eb7afU: // x = 0x1.dd6f5ep-7f 93 if (fputil::get_round() == FE_UPWARD) 94 return 0x1.d9fd86p-7f; 95 return 0x1.d9fd84p-7f; 96 case 0x41078febU: // x = 0x1.0f1fd6p+3f 97 if (fputil::get_round() != FE_UPWARD) 98 return 0x1.1fcbcep+1f; 99 break; 100 case 0x5cd69e88U: // x = 0x1.ad3d1p+58f 101 if (fputil::get_round() != FE_UPWARD) 102 return 0x1.45c146p+5f; 103 break; 104 case 0x65d890d3U: // x = 0x1.b121a6p+76f 105 if (fputil::get_round() == FE_TONEAREST) 106 return 0x1.a9a3f2p+5f; 107 break; 108 case 0x6f31a8ecU: // x = 0x1.6351d8p+95f 109 if (fputil::get_round() == FE_TONEAREST) 110 return 0x1.08b512p+6f; 111 break; 112 case 0x7a17f30aU: // x = 0x1.2fe614p+117f 113 if (fputil::get_round() != FE_UPWARD) 114 return 0x1.451436p+6f; 115 break; 116 case 0xbc4d092cU: // x = -0x1.9a1258p-7f 117 if (fputil::get_round() == FE_TONEAREST) 118 return -0x1.9ca8bep-7f; 119 break; 120 case 0xbc657728U: // x = -0x1.caee5p-7f 121 if (fputil::get_round() != FE_DOWNWARD) 122 return -0x1.ce2cccp-7f; 123 break; 124 case 0xbd1d20afU: // x = -0x1.3a415ep-5f 125 int round_mode = fputil::get_round(); 126 if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO) 127 return -0x1.40711p-5f; 128 return -0x1.407112p-5f; 129 } 130 131 return internal::log(xd + 1.0); 132 } 133 134 // Hard-to round cases. 135 switch (xbits.uintval()) { 136 case 0x35400003U: // x = 0x1.800006p-21f 137 if (fputil::get_round() == FE_TONEAREST) 138 return 0x1.7ffffep-21f; 139 break; 140 case 0x3710001bU: // x = 0x1.200036p-17f 141 if (fputil::get_round() == FE_TONEAREST) 142 return 0x1.1fffe6p-17f; 143 break; 144 case 0xb53ffffdU: // x = -0x1.7ffffap-21f 145 if (fputil::get_round() != FE_DOWNWARD) 146 return -0x1.800002p-21f; 147 break; 148 case 0xb70fffe5U: // x = -0x1.1fffcap-17f 149 if (fputil::get_round() != FE_DOWNWARD) 150 return -0x1.20001ap-17f; 151 break; 152 case 0xbb0ec8c4U: // x = -0x1.1d9188p-9f 153 if (fputil::get_round() == FE_TONEAREST) 154 return -0x1.1de14ap-9f; 155 break; 156 } 157 158 double r; 159 // Polymial generated with Sollya: 160 // > fpminimax(log(1 + x)/x, 5, [|D...|], [-2^-8; 2^-8]); 161 r = fputil::polyeval(xd, -0x1p-1, 0x1.5555555515551p-2, -0x1.ffffffff82bdap-3, 162 0x1.999b33348d3aep-3, -0x1.5556cae3adcc3p-3); 163 return static_cast<float>(fputil::multiply_add(r, xd * xd, xd)); 164 } 165 166 } // namespace __llvm_libc 167