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
log(double x)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