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