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