1 //===-- Double-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/log1p.h"
10 #include "src/__support/FPUtil/FEnvImpl.h"
11 #include "src/__support/FPUtil/FPBits.h"
12 #include "src/__support/FPUtil/PolyEval.h"
13 #include "src/__support/FPUtil/double_double.h"
14 #include "src/__support/FPUtil/dyadic_float.h"
15 #include "src/__support/FPUtil/multiply_add.h"
16 #include "src/__support/common.h"
17 #include "src/__support/integer_literals.h"
18 #include "src/__support/macros/config.h"
19 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
20
21 #include "common_constants.h"
22
23 namespace LIBC_NAMESPACE_DECL {
24
25 // 128-bit precision dyadic floating point numbers.
26 using Float128 = typename fputil::DyadicFloat<128>;
27
28 using LIBC_NAMESPACE::operator""_u128;
29
30 namespace {
31
32 // Extra errors from P is from using x^2 to reduce evaluation latency and
33 // directional rounding.
34 constexpr double P_ERR = 0x1.0p-49;
35
36 // log(2) with 128-bit precision generated by SageMath with:
37 // def format_hex(value):
38 // l = hex(value)[2:]
39 // n = 8
40 // x = [l[i:i + n] for i in range(0, len(l), n)]
41 // return "0x" + "'".join(x) + "_u128"
42 // (s, m, e) = RealField(128)(2).log().sign_mantissa_exponent();
43 // print(format_hex(m));
44 constexpr Float128 LOG_2(Sign::POS, /*exponent=*/-128, /*mantissa=*/
45 0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128);
46
47 // R1[i] = 2^-8 * nearestint( 2^8 / (1 + i * 2^-7) )
48 constexpr double R1[129] = {
49 0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.eap-1,
50 0x1.e6p-1, 0x1.e2p-1, 0x1.dep-1, 0x1.dap-1, 0x1.d8p-1, 0x1.d4p-1, 0x1.dp-1,
51 0x1.cep-1, 0x1.cap-1, 0x1.c8p-1, 0x1.c4p-1, 0x1.cp-1, 0x1.bep-1, 0x1.bap-1,
52 0x1.b8p-1, 0x1.b4p-1, 0x1.b2p-1, 0x1.bp-1, 0x1.acp-1, 0x1.aap-1, 0x1.a6p-1,
53 0x1.a4p-1, 0x1.a2p-1, 0x1.9ep-1, 0x1.9cp-1, 0x1.9ap-1, 0x1.98p-1, 0x1.94p-1,
54 0x1.92p-1, 0x1.9p-1, 0x1.8ep-1, 0x1.8ap-1, 0x1.88p-1, 0x1.86p-1, 0x1.84p-1,
55 0x1.82p-1, 0x1.8p-1, 0x1.7ep-1, 0x1.7ap-1, 0x1.78p-1, 0x1.76p-1, 0x1.74p-1,
56 0x1.72p-1, 0x1.7p-1, 0x1.6ep-1, 0x1.6cp-1, 0x1.6ap-1, 0x1.68p-1, 0x1.66p-1,
57 0x1.64p-1, 0x1.62p-1, 0x1.6p-1, 0x1.5ep-1, 0x1.5cp-1, 0x1.5ap-1, 0x1.58p-1,
58 0x1.58p-1, 0x1.56p-1, 0x1.54p-1, 0x1.52p-1, 0x1.5p-1, 0x1.4ep-1, 0x1.4cp-1,
59 0x1.4ap-1, 0x1.4ap-1, 0x1.48p-1, 0x1.46p-1, 0x1.44p-1, 0x1.42p-1, 0x1.42p-1,
60 0x1.4p-1, 0x1.3ep-1, 0x1.3cp-1, 0x1.3cp-1, 0x1.3ap-1, 0x1.38p-1, 0x1.36p-1,
61 0x1.36p-1, 0x1.34p-1, 0x1.32p-1, 0x1.3p-1, 0x1.3p-1, 0x1.2ep-1, 0x1.2cp-1,
62 0x1.2cp-1, 0x1.2ap-1, 0x1.28p-1, 0x1.28p-1, 0x1.26p-1, 0x1.24p-1, 0x1.24p-1,
63 0x1.22p-1, 0x1.2p-1, 0x1.2p-1, 0x1.1ep-1, 0x1.1cp-1, 0x1.1cp-1, 0x1.1ap-1,
64 0x1.1ap-1, 0x1.18p-1, 0x1.16p-1, 0x1.16p-1, 0x1.14p-1, 0x1.14p-1, 0x1.12p-1,
65 0x1.12p-1, 0x1.1p-1, 0x1.0ep-1, 0x1.0ep-1, 0x1.0cp-1, 0x1.0cp-1, 0x1.0ap-1,
66 0x1.0ap-1, 0x1.08p-1, 0x1.08p-1, 0x1.06p-1, 0x1.06p-1, 0x1.04p-1, 0x1.04p-1,
67 0x1.02p-1, 0x1.02p-1, 0x1p-1,
68 };
69
70 // Extra constants for exact range reduction when FMA instructions are not
71 // available:
72 // r * c - 1 for r = 2^-8 * nearestint( 2^8 / (1 + i * 2^-7))
73 // and c = 1 + i * 2^-7
74 // with i = 0..128.
75 [[maybe_unused]] constexpr double RCM1[129] = {
76 0.0, -0x1p-14, -0x1p-12, -0x1.2p-11, -0x1p-10, -0x1.9p-10,
77 0x1.fp-10, 0x1.28p-10, 0x1p-12, -0x1.9p-11, -0x1.fp-10, 0x1.2p-10,
78 -0x1p-12, -0x1.cp-10, 0x1.1p-10, -0x1.5p-11, 0x1p-9, 0x1p-14,
79 -0x1p-9, 0x1.ap-12, -0x1.ep-10, 0x1.8p-12, -0x1.1p-9, -0x1p-15,
80 0x1p-9, -0x1.ap-11, 0x1.1p-10, -0x1.f8p-10, -0x1p-12, 0x1.68p-10,
81 -0x1.fp-10, -0x1.cp-12, 0x1p-10, 0x1.3p-9, -0x1.6p-10, -0x1.4p-13,
82 0x1p-10, 0x1.0cp-9, -0x1.08p-9, -0x1.2p-10, -0x1p-12, 0x1.2p-11,
83 0x1.5p-10, 0x1p-9, 0x1.5p-9, -0x1.1cp-9, -0x1.cp-10, -0x1.58p-10,
84 -0x1p-10, -0x1.7p-11, -0x1p-11, -0x1.6p-12, -0x1p-12, -0x1.cp-13,
85 -0x1p-12, -0x1.6p-12, -0x1p-11, -0x1.7p-11, -0x1p-10, -0x1.58p-10,
86 -0x1.cp-10, -0x1.1cp-9, -0x1.6p-9, 0x1.5p-9, 0x1p-9, 0x1.5p-10,
87 0x1.2p-11, -0x1p-12, -0x1.2p-10, -0x1.08p-9, -0x1.88p-9, 0x1.0cp-9,
88 0x1p-10, -0x1.4p-13, -0x1.6p-10, -0x1.54p-9, 0x1.3p-9, 0x1p-10,
89 -0x1.cp-12, -0x1.fp-10, 0x1.8p-9, 0x1.68p-10, -0x1p-12, -0x1.f8p-10,
90 0x1.7p-9, 0x1.1p-10, -0x1.ap-11, -0x1.6p-9, 0x1p-9, -0x1p-15,
91 -0x1.1p-9, 0x1.48p-9, 0x1.8p-12, -0x1.ep-10, 0x1.6p-9, 0x1.ap-12,
92 -0x1p-9, 0x1.48p-9, 0x1p-14, -0x1.4p-9, 0x1p-9, -0x1.5p-11,
93 -0x1.bp-9, 0x1.1p-10, -0x1.cp-10, 0x1.54p-9, -0x1p-12, -0x1.9cp-9,
94 0x1.2p-10, -0x1.fp-10, 0x1.3p-9, -0x1.9p-11, 0x1.cp-9, 0x1p-12,
95 -0x1.88p-9, 0x1.28p-10, -0x1.2p-9, 0x1.fp-10, -0x1.9p-10, 0x1.4cp-9,
96 -0x1p-10, 0x1.9p-9, -0x1.2p-11, 0x1.c4p-9, -0x1p-12, 0x1.e8p-9,
97 -0x1p-14, 0x1.fcp-9, 0.0,
98 };
99
100 // Generated by Sollya with:
101 // for i from 0 to 128 do {
102 // r = 2^-8 * nearestint( 2^8 / (1 + i*2^-7) );
103 // b = nearestint(log(r)*2^43) * 2^-43;
104 // c = round(log(r) - b, D, RN);
105 // print("{", -c, ",", -b, "},");
106 // };
107 // We replace LOG_R1_DD[128] with log(1.0) == 0.0
108 constexpr fputil::DoubleDouble LOG_R1_DD[129] = {
109 {0.0, 0.0},
110 {-0x1.0c76b999d2be8p-46, 0x1.010157589p-7},
111 {-0x1.3dc5b06e2f7d2p-45, 0x1.0205658938p-6},
112 {-0x1.aa0ba325a0c34p-45, 0x1.8492528c9p-6},
113 {0x1.111c05cf1d753p-47, 0x1.0415d89e74p-5},
114 {-0x1.c167375bdfd28p-45, 0x1.466aed42ep-5},
115 {-0x1.29efbec19afa2p-47, 0x1.67c94f2d4cp-5},
116 {0x1.0fc1a353bb42ep-45, 0x1.aaef2d0fbp-5},
117 {-0x1.e113e4fc93b7bp-47, 0x1.eea31c006cp-5},
118 {-0x1.5325d560d9e9bp-45, 0x1.1973bd1466p-4},
119 {0x1.cc85ea5db4ed7p-45, 0x1.3bdf5a7d1ep-4},
120 {-0x1.53a2582f4e1efp-48, 0x1.4d3115d208p-4},
121 {0x1.c1e8da99ded32p-49, 0x1.700d30aeacp-4},
122 {0x1.3115c3abd47dap-45, 0x1.9335e5d594p-4},
123 {-0x1.e42b6b94407c8p-47, 0x1.a4e7640b1cp-4},
124 {0x1.646d1c65aacd3p-45, 0x1.c885801bc4p-4},
125 {0x1.a89401fa71733p-46, 0x1.da72763844p-4},
126 {-0x1.534d64fa10afdp-45, 0x1.fe89139dbep-4},
127 {0x1.1ef78ce2d07f2p-45, 0x1.1178e8227ep-3},
128 {0x1.ca78e44389934p-45, 0x1.1aa2b7e23fp-3},
129 {0x1.39d6ccb81b4a1p-47, 0x1.2d1610c868p-3},
130 {0x1.62fa8234b7289p-51, 0x1.365fcb0159p-3},
131 {0x1.5837954fdb678p-45, 0x1.4913d8333bp-3},
132 {0x1.633e8e5697dc7p-45, 0x1.527e5e4a1bp-3},
133 {-0x1.27023eb68981cp-46, 0x1.5bf406b544p-3},
134 {-0x1.5118de59c21e1p-45, 0x1.6f0128b757p-3},
135 {-0x1.c661070914305p-46, 0x1.7898d85445p-3},
136 {-0x1.73d54aae92cd1p-47, 0x1.8beafeb39p-3},
137 {0x1.7f22858a0ff6fp-47, 0x1.95a5adcf7p-3},
138 {0x1.9904d6865817ap-45, 0x1.9f6c407089p-3},
139 {-0x1.c358d4eace1aap-47, 0x1.b31d8575bdp-3},
140 {-0x1.d4bc4595412b6p-45, 0x1.bd087383bep-3},
141 {-0x1.1ec72c5962bd2p-48, 0x1.c6ffbc6f01p-3},
142 {-0x1.84a7e75b6f6e4p-47, 0x1.d1037f2656p-3},
143 {0x1.212276041f43p-51, 0x1.e530effe71p-3},
144 {-0x1.a211565bb8e11p-51, 0x1.ef5ade4ddp-3},
145 {0x1.bcbecca0cdf3p-46, 0x1.f991c6cb3bp-3},
146 {-0x1.6f08c1485e94ap-46, 0x1.01eae5626c8p-2},
147 {0x1.7188b163ceae9p-45, 0x1.0c42d67616p-2},
148 {-0x1.c210e63a5f01cp-45, 0x1.1178e8227e8p-2},
149 {0x1.b9acdf7a51681p-45, 0x1.16b5ccbacf8p-2},
150 {0x1.ca6ed5147bdb7p-45, 0x1.1bf99635a68p-2},
151 {0x1.a87deba46baeap-47, 0x1.214456d0eb8p-2},
152 {0x1.c93c1df5bb3b6p-45, 0x1.269621134d8p-2},
153 {0x1.a9cfa4a5004f4p-45, 0x1.2bef07cdc9p-2},
154 {0x1.16ecdb0f177c8p-46, 0x1.36b6776be1p-2},
155 {0x1.83b54b606bd5cp-46, 0x1.3c25277333p-2},
156 {0x1.8e436ec90e09dp-47, 0x1.419b423d5e8p-2},
157 {-0x1.f27ce0967d675p-45, 0x1.4718dc271c8p-2},
158 {-0x1.e20891b0ad8a4p-45, 0x1.4c9e09e173p-2},
159 {0x1.ebe708164c759p-45, 0x1.522ae0738ap-2},
160 {0x1.fadedee5d40efp-46, 0x1.57bf753c8dp-2},
161 {-0x1.a0b2a08a465dcp-47, 0x1.5d5bddf596p-2},
162 {-0x1.db623e731aep-45, 0x1.630030b3abp-2},
163 {0x1.0a0d32756ebap-45, 0x1.68ac83e9c68p-2},
164 {0x1.721657c222d87p-46, 0x1.6e60ee6af18p-2},
165 {0x1.d8b0949dc60b3p-45, 0x1.741d876c678p-2},
166 {0x1.9ec7d2efd1778p-45, 0x1.79e26687cf8p-2},
167 {-0x1.72090c812566ap-45, 0x1.7fafa3bd818p-2},
168 {0x1.fd56f3333778ap-45, 0x1.85855776dc8p-2},
169 {-0x1.05ae1e5e7047p-45, 0x1.8b639a88b3p-2},
170 {-0x1.766b52ee6307dp-46, 0x1.914a8635bf8p-2},
171 {-0x1.52313a502d9fp-46, 0x1.973a3431358p-2},
172 {-0x1.52313a502d9fp-46, 0x1.973a3431358p-2},
173 {-0x1.6279e10d0c0bp-45, 0x1.9d32bea15fp-2},
174 {0x1.3c6457f9d79f5p-45, 0x1.a33440224f8p-2},
175 {0x1.e36f2bea77a5dp-46, 0x1.a93ed3c8ad8p-2},
176 {-0x1.17cc552774458p-45, 0x1.af5295248dp-2},
177 {0x1.095252d841995p-46, 0x1.b56fa044628p-2},
178 {0x1.7d85bf40a666dp-45, 0x1.bb9611b80ep-2},
179 {0x1.cec807fe8e18p-45, 0x1.c1c60693fap-2},
180 {0x1.cec807fe8e18p-45, 0x1.c1c60693fap-2},
181 {-0x1.9b6ddc15249aep-45, 0x1.c7ff9c74558p-2},
182 {-0x1.797c33ec7a6bp-47, 0x1.ce42f180648p-2},
183 {0x1.35bafe9a767a8p-45, 0x1.d490246def8p-2},
184 {-0x1.ea42d60dc616ap-46, 0x1.dae75484c98p-2},
185 {-0x1.ea42d60dc616ap-46, 0x1.dae75484c98p-2},
186 {-0x1.326b207322938p-46, 0x1.e148a1a2728p-2},
187 {-0x1.465505372bd08p-45, 0x1.e7b42c3ddbp-2},
188 {0x1.f27f45a470251p-45, 0x1.ee2a156b41p-2},
189 {0x1.f27f45a470251p-45, 0x1.ee2a156b41p-2},
190 {0x1.2cde56f014a8bp-46, 0x1.f4aa7ee0318p-2},
191 {0x1.085fa3c164935p-47, 0x1.fb358af7a48p-2},
192 {-0x1.53ba3b1727b1cp-47, 0x1.00e5ae5b208p-1},
193 {-0x1.53ba3b1727b1cp-47, 0x1.00e5ae5b208p-1},
194 {-0x1.4c45fe79539ep-47, 0x1.04360be7604p-1},
195 {0x1.6812241edf5fdp-45, 0x1.078bf0533c4p-1},
196 {0x1.f486b887e7e27p-46, 0x1.0ae76e2d054p-1},
197 {0x1.f486b887e7e27p-46, 0x1.0ae76e2d054p-1},
198 {0x1.c299807801742p-46, 0x1.0e4898611ccp-1},
199 {-0x1.58647bb9ddcb2p-45, 0x1.11af823c75cp-1},
200 {-0x1.58647bb9ddcb2p-45, 0x1.11af823c75cp-1},
201 {-0x1.edd97a293ae49p-45, 0x1.151c3f6f298p-1},
202 {0x1.4cc4ef8ab465p-46, 0x1.188ee40f23cp-1},
203 {0x1.4cc4ef8ab465p-46, 0x1.188ee40f23cp-1},
204 {0x1.cacdeed70e667p-51, 0x1.1c07849ae6p-1},
205 {-0x1.a7242c9fe81d3p-45, 0x1.1f8635fc618p-1},
206 {-0x1.a7242c9fe81d3p-45, 0x1.1f8635fc618p-1},
207 {0x1.2fc066e48667bp-46, 0x1.230b0d8bebcp-1},
208 {-0x1.b61f10522625p-47, 0x1.269621134dcp-1},
209 {-0x1.b61f10522625p-47, 0x1.269621134dcp-1},
210 {0x1.06d2be797882dp-45, 0x1.2a2786d0ecp-1},
211 {-0x1.7a6e507b9dc11p-46, 0x1.2dbf557b0ep-1},
212 {-0x1.7a6e507b9dc11p-46, 0x1.2dbf557b0ep-1},
213 {-0x1.74e93c5a0ed9cp-45, 0x1.315da443408p-1},
214 {-0x1.74e93c5a0ed9cp-45, 0x1.315da443408p-1},
215 {0x1.0b83f9527e6acp-46, 0x1.35028ad9d8cp-1},
216 {-0x1.18b7abb5569a4p-45, 0x1.38ae2171978p-1},
217 {-0x1.18b7abb5569a4p-45, 0x1.38ae2171978p-1},
218 {-0x1.2b7367cfe13c2p-47, 0x1.3c6080c36cp-1},
219 {-0x1.2b7367cfe13c2p-47, 0x1.3c6080c36cp-1},
220 {-0x1.6ce7930f0c74cp-45, 0x1.4019c2125ccp-1},
221 {-0x1.6ce7930f0c74cp-45, 0x1.4019c2125ccp-1},
222 {-0x1.d984f481051f7p-48, 0x1.43d9ff2f924p-1},
223 {-0x1.2cb6af94d60aap-45, 0x1.47a1527e8a4p-1},
224 {-0x1.2cb6af94d60aap-45, 0x1.47a1527e8a4p-1},
225 {0x1.f7115ed4c541cp-49, 0x1.4b6fd6f970cp-1},
226 {0x1.f7115ed4c541cp-49, 0x1.4b6fd6f970cp-1},
227 {-0x1.e6c516d93b8fbp-45, 0x1.4f45a835a5p-1},
228 {-0x1.e6c516d93b8fbp-45, 0x1.4f45a835a5p-1},
229 {0x1.5ccc45d257531p-47, 0x1.5322e268678p-1},
230 {0x1.5ccc45d257531p-47, 0x1.5322e268678p-1},
231 {0x1.9980bff3303ddp-47, 0x1.5707a26bb8cp-1},
232 {0x1.9980bff3303ddp-47, 0x1.5707a26bb8cp-1},
233 {0x1.dfa63ac10c9fbp-45, 0x1.5af405c3648p-1},
234 {0x1.dfa63ac10c9fbp-45, 0x1.5af405c3648p-1},
235 {0x1.202380cda46bep-45, 0x1.5ee82aa2418p-1},
236 {0x1.202380cda46bep-45, 0x1.5ee82aa2418p-1},
237 {0.0, 0.0},
238 };
239
240 // Degree-7 minimax polynomial log(1 + v) ~ v - v^2 / 2 + ...
241 // generated by Sollya with:
242 // > P = fpminimax(log(1 + x)/x, 6, [|1, 1, D...|],
243 // [-0x1.69000000000edp-8, 0x1.7f00000000081p-8]);
244 constexpr double P_COEFFS[6] = {-0x1p-1,
245 0x1.5555555555166p-2,
246 -0x1.fffffffdb7746p-3,
247 0x1.99999a8718a6p-3,
248 -0x1.555874ce8ce22p-3,
249 0x1.24335555ddbe5p-3};
250
251 // -log(r1) with 128-bit precision generated by SageMath with:
252 //
253 // for i in range(129):
254 // r = 2^-8 * round( 2^8 / (1 + i*2^(-7)) );
255 // s, m, e = RealField(128)(r).log().sign_mantissa_exponent();
256 // print("{Sign::POS,", e, ", format_hex(m), "},");
257 constexpr Float128 LOG_R1[129] = {
258 {Sign::POS, 0, 0_u128},
259 {Sign::POS, -134, 0x8080abac'46f38946'662d417c'ed007a46_u128},
260 {Sign::POS, -133, 0x8102b2c4'9ac23a4f'91d082dc'e3ddcd38_u128},
261 {Sign::POS, -133, 0xc2492946'4655f45c'da5f3cc0'b3251dbd_u128},
262 {Sign::POS, -132, 0x820aec4f'3a222380'b9e3aea6'c444ef07_u128},
263 {Sign::POS, -132, 0xa33576a1'6f1f4c64'521016bd'904dc968_u128},
264 {Sign::POS, -132, 0xb3e4a796'a5dac208'27cca0bc'c06c2f92_u128},
265 {Sign::POS, -132, 0xd5779687'd887e0d1'a9dda170'56e45ed5_u128},
266 {Sign::POS, -132, 0xf7518e00'35c3dd83'606d8909'3278a939_u128},
267 {Sign::POS, -131, 0x8cb9de8a'32ab368a'a7c98595'30a45153_u128},
268 {Sign::POS, -131, 0x9defad3e'8f73217a'976d3b5b'45f6ca0b_u128},
269 {Sign::POS, -131, 0xa6988ae9'03f562ed'3e858f08'597b3a69_u128},
270 {Sign::POS, -131, 0xb8069857'560707a3'6a677b4c'8bec22e1_u128},
271 {Sign::POS, -131, 0xc99af2ea'ca4c4570'eaf51f66'692844ba_u128},
272 {Sign::POS, -131, 0xd273b205'8de1bd49'46bbf837'b4d320c6_u128},
273 {Sign::POS, -131, 0xe442c00d'e2591b47'196ab34c'e0bccd12_u128},
274 {Sign::POS, -131, 0xed393b1c'22351280'3f4e2e66'0317d55f_u128},
275 {Sign::POS, -131, 0xff4489ce'deab2ca6'c17bd40d'8d9291ec_u128},
276 {Sign::POS, -130, 0x88bc7411'3f23def1'9c5a0fe3'96f40f1e_u128},
277 {Sign::POS, -130, 0x8d515bf1'1fb94f1c'88713268'840cbcc0_u128},
278 {Sign::POS, -130, 0x968b0864'3409ceb6'65c0da50'6a088484_u128},
279 {Sign::POS, -130, 0x9b2fe580'ac80b17d'411a5b94'4aca8708_u128},
280 {Sign::POS, -130, 0xa489ec19'9dab06f2'a9fb6cf0'ecb411b7_u128},
281 {Sign::POS, -130, 0xa93f2f25'0dac67d1'cad2fb8d'48054ae0_u128},
282 {Sign::POS, -130, 0xadfa035a'a1ed8fdc'149767e4'10316d2c_u128},
283 {Sign::POS, -130, 0xb780945b'ab55dce4'34c7bc3d'32750fde_u128},
284 {Sign::POS, -130, 0xbc4c6c2a'226399ef'8f6ebcfb'2016a439_u128},
285 {Sign::POS, -130, 0xc5f57f59'c7f46155'aa8b6997'a402bf30_u128},
286 {Sign::POS, -130, 0xcad2d6e7'b80bf914'2c507fb7'a3d0bf6a_u128},
287 {Sign::POS, -130, 0xcfb62038'44b3209a'd0cb02f3'3f79c16c_u128},
288 {Sign::POS, -130, 0xd98ec2ba'de71e539'58a98f2a'd65bee9b_u128},
289 {Sign::POS, -130, 0xde8439c1'dec56877'4d57da94'5b5d0aaa_u128},
290 {Sign::POS, -130, 0xe37fde37'807b84e3'4e9a750b'6b68781d_u128},
291 {Sign::POS, -130, 0xe881bf93'2af3dac0'c524848e'3443e040_u128},
292 {Sign::POS, -130, 0xf29877ff'38809091'3b020fa1'820c9492_u128},
293 {Sign::POS, -130, 0xf7ad6f26'e7ff2ef7'54d2238f'75f969b1_u128},
294 {Sign::POS, -130, 0xfcc8e365'9d9bcbec'ca0cdf30'1431b60f_u128},
295 {Sign::POS, -129, 0x80f572b1'363487b9'f5bd0b5b'3479d5f4_u128},
296 {Sign::POS, -129, 0x86216b3b'0b17188b'163ceae8'8f720f1e_u128},
297 {Sign::POS, -129, 0x88bc7411'3f23def1'9c5a0fe3'96f40f1e_u128},
298 {Sign::POS, -129, 0x8b5ae65d'67db9acd'f7a51681'26a58b9a_u128},
299 {Sign::POS, -129, 0x8dfccb1a'd35ca6ed'5147bdb6'ddcaf59c_u128},
300 {Sign::POS, -129, 0x90a22b68'75c6a1f7'ae91aeba'609c8877_u128},
301 {Sign::POS, -129, 0x934b1089'a6dc93c1'df5bb3b6'0554e152_u128},
302 {Sign::POS, -129, 0x95f783e6'e49a9cfa'4a5004f3'ef063313_u128},
303 {Sign::POS, -129, 0x9b5b3bb5'f088b766'd878bbe3'd392be25_u128},
304 {Sign::POS, -129, 0x9e1293b9'998c1daa'5b035eae'273a855f_u128},
305 {Sign::POS, -129, 0xa0cda11e'af46390d'bb243827'3918db7e_u128},
306 {Sign::POS, -129, 0xa38c6e13'8e20d831'f698298a'dddd7f32_u128},
307 {Sign::POS, -129, 0xa64f04f0'b961df76'e4f5275c'2d15c21f_u128},
308 {Sign::POS, -129, 0xa9157039'c51ebe70'8164c759'686a2209_u128},
309 {Sign::POS, -129, 0xabdfba9e'468fd6f6'f72ea077'49ce6bd3_u128},
310 {Sign::POS, -129, 0xaeadeefa'caf97d35'7dd6e688'ebb13b03_u128},
311 {Sign::POS, -129, 0xb1801859'd56249dc'18ce51ff'f99479cd_u128},
312 {Sign::POS, -129, 0xb45641f4'e350a0d3'2756eba0'0bc33978_u128},
313 {Sign::POS, -129, 0xb7307735'78cb90b2'be1116c3'466beb6d_u128},
314 {Sign::POS, -129, 0xba0ec3b6'33dd8b09'49dc60b2'b059a60b_u128},
315 {Sign::POS, -129, 0xbcf13343'e7d9ec7d'2efd1778'1bb3afec_u128},
316 {Sign::POS, -129, 0xbfd7d1de'c0a8df6f'37eda996'244bccb0_u128},
317 {Sign::POS, -129, 0xc2c2abbb'6e5fd56f'33337789'd592e296_u128},
318 {Sign::POS, -129, 0xc5b1cd44'596fa51e'1a18fb8f'9f9ef280_u128},
319 {Sign::POS, -129, 0xc8a5431a'dfb44ca5'688ce7c1'a75e341a_u128},
320 {Sign::POS, -129, 0xcb9d1a18'9ab56e76'2d7e9307'c70c0668_u128},
321 {Sign::POS, -129, 0xcb9d1a18'9ab56e76'2d7e9307'c70c0668_u128},
322 {Sign::POS, -129, 0xce995f50'af69d861'ef2f3f4f'861ad6a9_u128},
323 {Sign::POS, -129, 0xd19a2011'27d3c645'7f9d79f5'1dcc7301_u128},
324 {Sign::POS, -129, 0xd49f69e4'56cf1b79'5f53bd2e'406e66e7_u128},
325 {Sign::POS, -129, 0xd7a94a92'466e833a'ad88bba7'd0cee8e0_u128},
326 {Sign::POS, -129, 0xdab7d022'31484a92'96c20cca'6efe2ac5_u128},
327 {Sign::POS, -129, 0xddcb08dc'0717d85b'f40a666c'87842843_u128},
328 {Sign::POS, -129, 0xe0e30349'fd1cec80'7fe8e180'2aba24d6_u128},
329 {Sign::POS, -129, 0xe0e30349'fd1cec80'7fe8e180'2aba24d6_u128},
330 {Sign::POS, -129, 0xe3ffce3a'2aa64922'3eadb651'b49ac53a_u128},
331 {Sign::POS, -129, 0xe72178c0'323a1a0f'304e1653'e71d9973_u128},
332 {Sign::POS, -129, 0xea481236'f7d35baf'e9a767a8'0d6d97e8_u128},
333 {Sign::POS, -129, 0xed73aa42'64b0ade9'4f91cf4b'33e42998_u128},
334 {Sign::POS, -129, 0xed73aa42'64b0ade9'4f91cf4b'33e42998_u128},
335 {Sign::POS, -129, 0xf0a450d1'39366ca6'fc66eb64'08ff6433_u128},
336 {Sign::POS, -129, 0xf3da161e'ed6b9aaf'ac8d42f7'8d3e65d3_u128},
337 {Sign::POS, -129, 0xf7150ab5'a09f27f4'5a470250'd40ebe90_u128},
338 {Sign::POS, -129, 0xf7150ab5'a09f27f4'5a470250'd40ebe90_u128},
339 {Sign::POS, -129, 0xfa553f70'18c966f2'b780a545'a1b54dcf_u128},
340 {Sign::POS, -129, 0xfd9ac57b'd244217e'8f05924d'258c14c5_u128},
341 {Sign::POS, -128, 0x8072d72d'903d588b'89d1b09c'70c4010a_u128},
342 {Sign::POS, -128, 0x8072d72d'903d588b'89d1b09c'70c4010a_u128},
343 {Sign::POS, -128, 0x821b05f3'b01d6774'030d58c3'f7e2ea1f_u128},
344 {Sign::POS, -128, 0x83c5f829'9e2b4091'20f6fafe'8fbb68b9_u128},
345 {Sign::POS, -128, 0x8573b716'82a7d21a'e21f9f89'c1ab80b2_u128},
346 {Sign::POS, -128, 0x8573b716'82a7d21a'e21f9f89'c1ab80b2_u128},
347 {Sign::POS, -128, 0x87244c30'8e670a66'01e005d0'6dbfa8f8_u128},
348 {Sign::POS, -128, 0x88d7c11e'3ad53cdc'223111a7'07b6de2c_u128},
349 {Sign::POS, -128, 0x88d7c11e'3ad53cdc'223111a7'07b6de2c_u128},
350 {Sign::POS, -128, 0x8a8e1fb7'94b09134'2eb628db'a173c82d_u128},
351 {Sign::POS, -128, 0x8c477207'91e53313'be2ad194'15fe25a5_u128},
352 {Sign::POS, -128, 0x8c477207'91e53313'be2ad194'15fe25a5_u128},
353 {Sign::POS, -128, 0x8e03c24d'73003959'bddae1cc'ce247838_u128},
354 {Sign::POS, -128, 0x8fc31afe'30b2c6de'9b00bf16'7e95da67_u128},
355 {Sign::POS, -128, 0x8fc31afe'30b2c6de'9b00bf16'7e95da67_u128},
356 {Sign::POS, -128, 0x918586c5'f5e4bf01'9b92199e'd1a4bab1_u128},
357 {Sign::POS, -128, 0x934b1089'a6dc93c1'df5bb3b6'0554e152_u128},
358 {Sign::POS, -128, 0x934b1089'a6dc93c1'df5bb3b6'0554e152_u128},
359 {Sign::POS, -128, 0x9513c368'76083695'f3cbc416'a2418012_u128},
360 {Sign::POS, -128, 0x96dfaabd'86fa1646'be1188fb'c94e2f15_u128},
361 {Sign::POS, -128, 0x96dfaabd'86fa1646'be1188fb'c94e2f15_u128},
362 {Sign::POS, -128, 0x98aed221'a03458b6'1d2f8932'1647b358_u128},
363 {Sign::POS, -128, 0x98aed221'a03458b6'1d2f8932'1647b358_u128},
364 {Sign::POS, -128, 0x9a81456c'ec642e0f'e549f9aa'ea3cb5e1_u128},
365 {Sign::POS, -128, 0x9c5710b8'cbb73a42'a2554b2d'd4619e63_u128},
366 {Sign::POS, -128, 0x9c5710b8'cbb73a42'a2554b2d'd4619e63_u128},
367 {Sign::POS, -128, 0x9e304061'b5fda919'30603d87'b6df81ad_u128},
368 {Sign::POS, -128, 0x9e304061'b5fda919'30603d87'b6df81ad_u128},
369 {Sign::POS, -128, 0xa00ce109'2e5498c3'67879c5a'30cd1242_u128},
370 {Sign::POS, -128, 0xa00ce109'2e5498c3'67879c5a'30cd1242_u128},
371 {Sign::POS, -128, 0xa1ecff97'c91e267b'0b7efae0'8e597e16_u128},
372 {Sign::POS, -128, 0xa3d0a93f'45169a4a'83594fab'088c0d65_u128},
373 {Sign::POS, -128, 0xa3d0a93f'45169a4a'83594fab'088c0d65_u128},
374 {Sign::POS, -128, 0xa5b7eb7c'b860fb88'af6a62a0'dec6e073_u128},
375 {Sign::POS, -128, 0xa5b7eb7c'b860fb88'af6a62a0'dec6e073_u128},
376 {Sign::POS, -128, 0xa7a2d41a'd270c9d7'49362382'a768847a_u128},
377 {Sign::POS, -128, 0xa7a2d41a'd270c9d7'49362382'a768847a_u128},
378 {Sign::POS, -128, 0xa9917134'33c2b998'8ba4aea6'14d05701_u128},
379 {Sign::POS, -128, 0xa9917134'33c2b998'8ba4aea6'14d05701_u128},
380 {Sign::POS, -128, 0xab83d135'dc633301'7fe6607b'a902ef3c_u128},
381 {Sign::POS, -128, 0xab83d135'dc633301'7fe6607b'a902ef3c_u128},
382 {Sign::POS, -128, 0xad7a02e1'b24efd31'd60864fd'949b4bd3_u128},
383 {Sign::POS, -128, 0xad7a02e1'b24efd31'd60864fd'949b4bd3_u128},
384 {Sign::POS, -128, 0xaf741551'20c9011c'066d235e'e63073dd_u128},
385 {Sign::POS, -128, 0xaf741551'20c9011c'066d235e'e63073dd_u128},
386 {Sign::POS, 0, 0_u128},
387 };
388
389 // Logarithm range reduction - Step 2:
390 // s(k) = 2^-18 round( 2^18 / (1 + k*2^-14) ) - 1 for k = -91 .. 96
391 // Output range:
392 // [-0x1.1037c00000040271p-15 , 0x1.108480000008096cp-15]
393 constexpr double S2[198] = {
394 0x1.6ep-8, 0x1.6ap-8, 0x1.66p-8, 0x1.62p-8, 0x1.5dcp-8,
395 0x1.59cp-8, 0x1.55cp-8, 0x1.51cp-8, 0x1.4dcp-8, 0x1.49cp-8,
396 0x1.458p-8, 0x1.418p-8, 0x1.3d8p-8, 0x1.398p-8, 0x1.358p-8,
397 0x1.318p-8, 0x1.2d8p-8, 0x1.294p-8, 0x1.254p-8, 0x1.214p-8,
398 0x1.1d4p-8, 0x1.194p-8, 0x1.154p-8, 0x1.114p-8, 0x1.0dp-8,
399 0x1.09p-8, 0x1.05p-8, 0x1.01p-8, 0x1.fap-9, 0x1.f2p-9,
400 0x1.eap-9, 0x1.e2p-9, 0x1.d98p-9, 0x1.d18p-9, 0x1.c98p-9,
401 0x1.c18p-9, 0x1.b98p-9, 0x1.b18p-9, 0x1.a98p-9, 0x1.a18p-9,
402 0x1.998p-9, 0x1.91p-9, 0x1.89p-9, 0x1.81p-9, 0x1.79p-9,
403 0x1.71p-9, 0x1.69p-9, 0x1.61p-9, 0x1.59p-9, 0x1.51p-9,
404 0x1.49p-9, 0x1.41p-9, 0x1.388p-9, 0x1.308p-9, 0x1.288p-9,
405 0x1.208p-9, 0x1.188p-9, 0x1.108p-9, 0x1.088p-9, 0x1.008p-9,
406 0x1.f1p-10, 0x1.e1p-10, 0x1.d1p-10, 0x1.c1p-10, 0x1.b1p-10,
407 0x1.a1p-10, 0x1.91p-10, 0x1.81p-10, 0x1.71p-10, 0x1.6p-10,
408 0x1.5p-10, 0x1.4p-10, 0x1.3p-10, 0x1.2p-10, 0x1.1p-10,
409 0x1p-10, 0x1.ep-11, 0x1.cp-11, 0x1.ap-11, 0x1.8p-11,
410 0x1.6p-11, 0x1.4p-11, 0x1.2p-11, 0x1p-11, 0x1.cp-12,
411 0x1.8p-12, 0x1.4p-12, 0x1p-12, 0x1.8p-13, 0x1p-13,
412 0x1p-14, 0.0, -0x1p-14, -0x1p-13, -0x1.8p-13,
413 -0x1p-12, -0x1.4p-12, -0x1.8p-12, -0x1.cp-12, -0x1p-11,
414 -0x1.2p-11, -0x1.4p-11, -0x1.6p-11, -0x1.8p-11, -0x1.ap-11,
415 -0x1.cp-11, -0x1.ep-11, -0x1p-10, -0x1.1p-10, -0x1.2p-10,
416 -0x1.3p-10, -0x1.4p-10, -0x1.5p-10, -0x1.6p-10, -0x1.6fp-10,
417 -0x1.7fp-10, -0x1.8fp-10, -0x1.9fp-10, -0x1.afp-10, -0x1.bfp-10,
418 -0x1.cfp-10, -0x1.dfp-10, -0x1.efp-10, -0x1.ffp-10, -0x1.078p-9,
419 -0x1.0f8p-9, -0x1.178p-9, -0x1.1f8p-9, -0x1.278p-9, -0x1.2f8p-9,
420 -0x1.378p-9, -0x1.3fp-9, -0x1.47p-9, -0x1.4fp-9, -0x1.57p-9,
421 -0x1.5fp-9, -0x1.67p-9, -0x1.6fp-9, -0x1.77p-9, -0x1.7fp-9,
422 -0x1.87p-9, -0x1.8fp-9, -0x1.968p-9, -0x1.9e8p-9, -0x1.a68p-9,
423 -0x1.ae8p-9, -0x1.b68p-9, -0x1.be8p-9, -0x1.c68p-9, -0x1.ce8p-9,
424 -0x1.d68p-9, -0x1.dep-9, -0x1.e6p-9, -0x1.eep-9, -0x1.f6p-9,
425 -0x1.fep-9, -0x1.03p-8, -0x1.07p-8, -0x1.0bp-8, -0x1.0fp-8,
426 -0x1.12cp-8, -0x1.16cp-8, -0x1.1acp-8, -0x1.1ecp-8, -0x1.22cp-8,
427 -0x1.26cp-8, -0x1.2acp-8, -0x1.2e8p-8, -0x1.328p-8, -0x1.368p-8,
428 -0x1.3a8p-8, -0x1.3e8p-8, -0x1.428p-8, -0x1.464p-8, -0x1.4a4p-8,
429 -0x1.4e4p-8, -0x1.524p-8, -0x1.564p-8, -0x1.5a4p-8, -0x1.5ep-8,
430 -0x1.62p-8, -0x1.66p-8, -0x1.6ap-8, -0x1.6ep-8, -0x1.72p-8,
431 -0x1.75cp-8, -0x1.79cp-8, -0x1.7dcp-8,
432 };
433
434 // -log(r) for the second step, generated by SageMath with:
435 //
436 // for i in range(-91, 97):
437 // r = 2^-18 * round( 2^18 / (1 + i*2^(-14)) );
438 // s, m, e = RealField(128)(r).log().sign_mantissa_exponent();
439 // print("{Sign::POS," if (s == -1) else "{Sign::NEG,", e, ",
440 // format_hex(m), "},");
441 constexpr Float128 LOG_R2[198] = {
442 {Sign::NEG, -135, 0xb67dab2a'1a5742a4'a0e061c5'f7431c5e_u128},
443 {Sign::NEG, -135, 0xb4807f24'af682939'5d5bfe7b'969ed6ec_u128},
444 {Sign::NEG, -135, 0xb2834b35'b4d54d5f'4d08702d'dfabc23f_u128},
445 {Sign::NEG, -135, 0xb0860f5c'eba9be95'd4d36650'8b9953df_u128},
446 {Sign::NEG, -135, 0xae68f71a'a09e8847'ac18a289'f8f214a9_u128},
447 {Sign::NEG, -135, 0xac6baaee'd676e8f1'd5b42054'abb88c45_u128},
448 {Sign::NEG, -135, 0xaa6e56d8'7cd632d6'09809d58'ee484964_u128},
449 {Sign::NEG, -135, 0xa870fad7'54bb8791'b9e6fc7c'72f06d73_u128},
450 {Sign::NEG, -135, 0xa67396eb'1f231892'6f78d6d0'105c00e2_u128},
451 {Sign::NEG, -135, 0xa4762b13'9d0626e7'028f7126'29209148_u128},
452 {Sign::NEG, -135, 0xa258dfd1'0aedaa67'c98d898e'f172df02_u128},
453 {Sign::NEG, -135, 0xa05b63a3'73e60a83'fcc37c3c'3062bfa1_u128},
454 {Sign::NEG, -135, 0x9e5ddf89'cf42f501'3eb450db'05763c36_u128},
455 {Sign::NEG, -135, 0x9c605383'ddf1b88c'7146a86f'd458b775_u128},
456 {Sign::NEG, -135, 0x9a62bf91'60dcb286'c20a0c92'81474436_u128},
457 {Sign::NEG, -135, 0x986523b2'18eb4ed6'cdc57316'ec4aebc3_u128},
458 {Sign::NEG, -135, 0x96677fe5'c70207b9'c060dad7'4cef4273_u128},
459 {Sign::NEG, -135, 0x9449f92d'2ff44633'ed8def1a'3e433499_u128},
460 {Sign::NEG, -135, 0x924c4507'3220b5e0'3ce7a1f8'5c27b4fc_u128},
461 {Sign::NEG, -135, 0x904e88f3'68fea63f'f2ca8934'49f7f2cb_u128},
462 {Sign::NEG, -135, 0x8e50c4f1'956699ed'8d77d9fa'bd2853cf_u128},
463 {Sign::NEG, -135, 0x8c52f901'782e20ec'93e828d7'5b58ded4_u128},
464 {Sign::NEG, -135, 0x8a552522'd227d87a'9f9605b0'53c5acf0_u128},
465 {Sign::NEG, -135, 0x88574955'64236ae0'62a14939'3bca7241_u128},
466 {Sign::NEG, -135, 0x86398719'b66bac7c'aea6b56c'e89203d4_u128},
467 {Sign::NEG, -135, 0x843b9aef'044e4dcc'0242bd86'd00609b2_u128},
468 {Sign::NEG, -135, 0x823da6d4'c89c6927'daabf927'74bac84e_u128},
469 {Sign::NEG, -135, 0x803faaca'c419abf2'a1c6f3fc'242ef8d0_u128},
470 {Sign::NEG, -136, 0xfc834da1'6f0d9f57'a225ebc0'2e6d9dd4_u128},
471 {Sign::NEG, -136, 0xf88735cc'c7433381'c33f6ad3'40ae18a9_u128},
472 {Sign::NEG, -136, 0xf48b0e17'1249b6bc'70b2a4d3'8a242244_u128},
473 {Sign::NEG, -136, 0xf08ed67f'd190e280'1d548190'48b811b0_u128},
474 {Sign::NEG, -136, 0xec52ca07'ed95f236'9c21b650'afe9ede0_u128},
475 {Sign::NEG, -136, 0xe85671ad'ecd28aac'935519c9'6d30e463_u128},
476 {Sign::NEG, -136, 0xe45a0970'dc912ca7'ba88f6f2'e2672cfe_u128},
477 {Sign::NEG, -136, 0xe05d9150'3e298bc8'0b1a8b84'657ae069_u128},
478 {Sign::NEG, -136, 0xdc61094b'92ed70ef'ea3bff8d'197b20a1_u128},
479 {Sign::NEG, -136, 0xd8647162'5c28b9e5'cdbb931d'6fecc249_u128},
480 {Sign::NEG, -136, 0xd467c994'1b2158f5'd971d560'd5f00820_u128},
481 {Sign::NEG, -136, 0xd06b11e0'51175493'75563561'244c090b_u128},
482 {Sign::NEG, -136, 0xcc6e4a46'7f44c6fa'dc393c9a'3f3b380f_u128},
483 {Sign::NEG, -136, 0xc831a4c6'f6fa709d'e6abe6e9'e4ee2096_u128},
484 {Sign::NEG, -136, 0xc434bc61'24a0f16e'3ce3c822'8583a66e_u128},
485 {Sign::NEG, -136, 0xc037c413'c61bfd93'b96a79f5'c5a4963a_u128},
486 {Sign::NEG, -136, 0xbc3abbde'5c8d9bde'aaef2733'7008679f_u128},
487 {Sign::NEG, -136, 0xb83da3c0'6911e509'a49a3fca'ddc8bc5a_u128},
488 {Sign::NEG, -136, 0xb4407bb9'6cbf035a'e0254feb'785362fa_u128},
489 {Sign::NEG, -136, 0xb04343c8'e8a53245'9893a4e2'5ab9dc95_u128},
490 {Sign::NEG, -136, 0xac45fbee'5dcebe0b'5d8b0f40'a3708915_u128},
491 {Sign::NEG, -136, 0xa848a429'4d40035d'5f4c11c2'c7a58c69_u128},
492 {Sign::NEG, -136, 0xa44b3c79'37f76efd'b348cc5d'f706ffba_u128},
493 {Sign::NEG, -136, 0xa04dc4dd'9eed7d60'9159f2c5'5a18befd_u128},
494 {Sign::NEG, -136, 0x9c106456'3058bef3'bdfdee41'fe6a5a02_u128},
495 {Sign::NEG, -136, 0x9812cbe3'46475a24'4580ddf8'9853254d_u128},
496 {Sign::NEG, -136, 0x94152383'53489ffb'ac75e10d'61fc3ee8_u128},
497 {Sign::NEG, -136, 0x90176b35'd83ce8e2'cad9b30b'29736155_u128},
498 {Sign::NEG, -136, 0x8c19a2fa'55fe9b14'6f881deb'98fc45f3_u128},
499 {Sign::NEG, -136, 0x881bcad0'4d622a3e'70a04b63'b7248c96_u128},
500 {Sign::NEG, -136, 0x841de2b7'3f361722'b4823fb4'8035eddd_u128},
501 {Sign::NEG, -136, 0x801feaae'ac42ef38'3364ccb5'b13cd47f_u128},
502 {Sign::NEG, -137, 0xf843c56c'2a969897'e306977b'049f0ad5_u128},
503 {Sign::NEG, -137, 0xf0479599'f617a843'e3c4d9e9'619bc045_u128},
504 {Sign::NEG, -137, 0xe84b45e5'bc76702c'4356d525'b5e6432d_u128},
505 {Sign::NEG, -137, 0xe04ed64e'7f14697a'7839dcd7'989339ab_u128},
506 {Sign::NEG, -137, 0xd85246d3'3f47230b'4e21f045'ecb76f23_u128},
507 {Sign::NEG, -137, 0xd0559772'fe5840b0'902e248d'd4ba9b28_u128},
508 {Sign::NEG, -137, 0xc858c82c'bd857a72'a4444906'7ef92e01_u128},
509 {Sign::NEG, -137, 0xc05bd8ff'7e009bd2'17926207'cc22e4e6_u128},
510 {Sign::NEG, -137, 0xb85ec9ea'40ef8309'1c349622'f3fa5d82_u128},
511 {Sign::NEG, -137, 0xafe1c6ec'e1a058dd'97fa2fd0'c9dc723e_u128},
512 {Sign::NEG, -137, 0xa7e47606'048b1a65'983e8089'7cf1e60f_u128},
513 {Sign::NEG, -137, 0x9fe70534'1d236102'7199cd06'ae5d39b3_u128},
514 {Sign::NEG, -137, 0x97e97476'2c5e8f58'43cd18a7'2a051a96_u128},
515 {Sign::NEG, -137, 0x8febc3cb'332616ff'7b6d1248'c3e1fd40_u128},
516 {Sign::NEG, -137, 0x87edf332'325777c5'f5572a88'14c703af_u128},
517 {Sign::NEG, -138, 0xffe00554'55887de0'26828c92'649a3a39_u128},
518 {Sign::NEG, -138, 0xefe3e464'3a640cf3'82c550bd'1216d82a_u128},
519 {Sign::NEG, -138, 0xdfe78392'14b4e8ae'da6959f7'f0e01bf0_u128},
520 {Sign::NEG, -138, 0xcfeae2db'e5d6736d'da93e2fa'85a8f214_u128},
521 {Sign::NEG, -138, 0xbfee023f'af0c2480'b47505bf'a5a03b06_u128},
522 {Sign::NEG, -138, 0xaff0e1bb'718186ad'b1475a51'80a43520_u128},
523 {Sign::NEG, -138, 0x9ff3814d'2e4a36b2'a8740b91'c95df537_u128},
524 {Sign::NEG, -138, 0x8ff5e0f2'e661e1c6'57d895d3'5921b59c_u128},
525 {Sign::NEG, -139, 0xfff00155'35588833'3c56c598'c659c2a3_u128},
526 {Sign::NEG, -139, 0xdff3c0e4'97ea4eb1'2ef8ec33'ed9d782a_u128},
527 {Sign::NEG, -139, 0xbff7008f'f5e0c257'379eba7e'6465ff63_u128},
528 {Sign::NEG, -139, 0x9ff9c053'5073a370'3f972b78'3fcab757_u128},
529 {Sign::NEG, -140, 0xfff80055'51558885'de026e27'1ee0549d_u128},
530 {Sign::NEG, -140, 0xbffb8023'febc0c25'eceb47ea'01f6c632_u128},
531 {Sign::NEG, -141, 0xfffc0015'54d55888'7333c578'57e1ed52_u128},
532 {Sign::NEG, -142, 0xfffe0005'55455588'87dde026'fa704374_u128},
533 {Sign::POS, 0, 0_u128},
534 {Sign::POS, -141, 0x80010002'aab2aac4'44999abe'2fe2cc65_u128},
535 {Sign::POS, -140, 0x8002000a'aaeaac44'4eef3815'81464ccb_u128},
536 {Sign::POS, -140, 0xc0048024'01440c26'dfeb4850'85f6f454_u128},
537 {Sign::POS, -139, 0x8004002a'acaac445'99abe3be'3a1c6e93_u128},
538 {Sign::POS, -139, 0xa0064053'5a37a37a'6bc1e20e'ac8448b4_u128},
539 {Sign::POS, -139, 0xc0090090'0a20c275'979eedc0'64c242fd_u128},
540 {Sign::POS, -139, 0xe00c40e4'bd6e4efd'c72446cc'1bf728bd_u128},
541 {Sign::POS, -138, 0x800800aa'baac446e'f381b821'bbb569e5_u128},
542 {Sign::POS, -138, 0x900a20f3'19a3e273'569b26aa'a485ea5c_u128},
543 {Sign::POS, -138, 0xa00c814d'7c6a37f8'2dcf56c8'3c80b028_u128},
544 {Sign::POS, -138, 0xb00f21bb'e3e388ee'5f697682'84463b9b_u128},
545 {Sign::POS, -138, 0xc0120240'510c284c'b48ea6c0'5e2773a1_u128},
546 {Sign::POS, -138, 0xd01522dc'c4f87991'14d9d761'96d8043a_u128},
547 {Sign::POS, -138, 0xe0188393'40d4f241'e016a611'a4415d72_u128},
548 {Sign::POS, -138, 0xf01c2465'c5e61b6f'661e135f'49a47c40_u128},
549 {Sign::POS, -137, 0x801002ab'2ac4499a'be6bf0fa'435e8383_u128},
550 {Sign::POS, -137, 0x88121333'7898871e'9a31ba0c'bc030353_u128},
551 {Sign::POS, -137, 0x901443cc'cd362c9f'54b57dfe'0c4c840f_u128},
552 {Sign::POS, -137, 0x98169478'296fad41'7ad1e9c3'15328f7e_u128},
553 {Sign::POS, -137, 0xa0190536'8e2389b3'1f3f686c'f3d6be22_u128},
554 {Sign::POS, -137, 0xa81b9608'fc3c50ec'f105b66e'c4703ede_u128},
555 {Sign::POS, -137, 0xb01e46f0'74b0a0f3'610848c6'8df4d233_u128},
556 {Sign::POS, -137, 0xb7a0e9ed'7613acb0'2e0efddf'33a20464_u128},
557 {Sign::POS, -137, 0xbfa3d900'8e042ffb'c2cdb3c7'50f127b4_u128},
558 {Sign::POS, -137, 0xc7a6e82b'a36a7073'bd953378'6d3f4c49_u128},
559 {Sign::POS, -137, 0xcfaa176f'b76c8eb1'82e237c9'a4d450e3_u128},
560 {Sign::POS, -137, 0xd7ad66cd'cb3cbe14'c00b46a4'd0e3dfd0_u128},
561 {Sign::POS, -137, 0xdfb0d646'e0194584'ea999c0d'f8546710_u128},
562 {Sign::POS, -137, 0xe7b465db'f74c8032'cec6c2a9'ad974f4f_u128},
563 {Sign::POS, -137, 0xefb8158e'122cde5a'2d2045da'1570a07c_u128},
564 {Sign::POS, -137, 0xf7bbe55e'321ce603'6752e9b2'381e3edc_u128},
565 {Sign::POS, -137, 0xffbfd54d'588b33c5'3c1ed527'28e00e40_u128},
566 {Sign::POS, -136, 0x83e1f2ae'43793dc3'493b0d87'3fb9a340_u128},
567 {Sign::POS, -136, 0x87e40ac6'5f6cc4a0'29e38750'c9d26893_u128},
568 {Sign::POS, -136, 0x8be632ef'80e9a0df'aab9e832'7258ac3f_u128},
569 {Sign::POS, -136, 0x8fe86b2a'28bf51b3'28bc403d'8a5f3c63_u128},
570 {Sign::POS, -136, 0x93eab376'd7c36377'f720c1c9'7227fcdc_u128},
571 {Sign::POS, -136, 0x97ed0bd6'0ed17018'6ad9a3e3'd11b66c1_u128},
572 {Sign::POS, -136, 0x9bef7448'4ecb1f6c'edb27b79'c90b4019_u128},
573 {Sign::POS, -136, 0x9fb1c4cd'27012e19'a092a0d7'ab21722a_u128},
574 {Sign::POS, -136, 0xa3b44c65'b71c2d85'535d52f0'939a4d02_u128},
575 {Sign::POS, -136, 0xa7b6e412'cadcb3dc'90a57e11'edc1864e_u128},
576 {Sign::POS, -136, 0xabb98bd4'e33c4381'68e9c901'60031159_u128},
577 {Sign::POS, -136, 0xafbc43ac'813a6ea3'bf60594f'929adeb8_u128},
578 {Sign::POS, -136, 0xb3bf0b9a'25dcd7a2'8a421588'86775205_u128},
579 {Sign::POS, -136, 0xb7c1e39e'522f316d'1ab45417'663dee9e_u128},
580 {Sign::POS, -136, 0xbbc4cbb9'87433fe4'6c51ae3c'e1aea68a_u128},
581 {Sign::POS, -136, 0xbfc7c3ec'4630d83c'7c52ae8b'40ebabb7_u128},
582 {Sign::POS, -136, 0xc3cacc37'1015e15d'a857126f'7cfaaa67_u128},
583 {Sign::POS, -136, 0xc7cde49a'66165446'14d05662'cd29464a_u128},
584 {Sign::POS, -136, 0xcb90da16'44d29bb7'8379db06'ef3cd6bb_u128},
585 {Sign::POS, -136, 0xcf9411aa'99ddb7de'9025f4c6'7dd38bb6_u128},
586 {Sign::POS, -136, 0xd3975958'f681086d'd6f8a61c'892032ee_u128},
587 {Sign::POS, -136, 0xd79ab121'dbf8714c'9a2f20b4'e2332d47_u128},
588 {Sign::POS, -136, 0xdb9e1905'cb85ea59'3c767d61'f51d375b_u128},
589 {Sign::POS, -136, 0xdfa19105'46717fca'd4b2bd65'bb25493c_u128},
590 {Sign::POS, -136, 0xe3a51920'ce095292'c96c1254'a30ef91f_u128},
591 {Sign::POS, -136, 0xe7a8b158'e3a198be'73e324ce'0946b214_u128},
592 {Sign::POS, -136, 0xebac59ae'08949dd8'cacd125a'12bac62c_u128},
593 {Sign::POS, -136, 0xef6fd620'b2b7a503'cafdc272'27b71eaa_u128},
594 {Sign::POS, -136, 0xf3739daf'959aaafc'688d4282'f6026aa3_u128},
595 {Sign::POS, -136, 0xf777755d'03f4e0b6'e54e9e38'04464cdd_u128},
596 {Sign::POS, -136, 0xfb7b5d29'7f388a12'cb78b383'f4b59dce_u128},
597 {Sign::POS, -136, 0xff7f5515'88de024f'ee055fc5'15062c04_u128},
598 {Sign::POS, -135, 0x81c1ae90'd131de38'207812b4'3382acdd_u128},
599 {Sign::POS, -135, 0x83c3baa7'26a721cc'dc90c4c4'b61f3a87_u128},
600 {Sign::POS, -135, 0x85c5cece'05941dbc'1a03f13f'b2c978b1_u128},
601 {Sign::POS, -135, 0x87c7eb05'aec1304f'b36f282e'83a7dc36_u128},
602 {Sign::POS, -135, 0x89a9eccd'56a980c0'd82a4661'6d4c393f_u128},
603 {Sign::POS, -135, 0x8bac18a6'40185360'bc6ff847'13c9babd_u128},
604 {Sign::POS, -135, 0x8dae4c90'b22574f4'9f7942a5'16fc2d8a_u128},
605 {Sign::POS, -135, 0x8fb0888c'eda546ab'15e50cfd'9b29b427_u128},
606 {Sign::POS, -135, 0x91b2cc9b'336f3718'9f465296'ae7dd49a_u128},
607 {Sign::POS, -135, 0x93b518bb'c45dc268'b49c1eb9'b348e6e4_u128},
608 {Sign::POS, -135, 0x95b76cee'e14e728e'daa320cd'64c9d9c7_u128},
609 {Sign::POS, -135, 0x9799a333'de49b963'75a91950'ffe1e3b5_u128},
610 {Sign::POS, -135, 0x999c070b'a32068cd'5c6abcbf'43f03f14_u128},
611 {Sign::POS, -135, 0x9b9e72f6'b295ad4f'5a9e7f26'5d1ed157_u128},
612 {Sign::POS, -135, 0x9da0e6f5'4d9318fd'efeb98d0'2a195c17_u128},
613 {Sign::POS, -135, 0x9fa36307'b5054ca8'2aa503a3'110ab5a7_u128},
614 {Sign::POS, -135, 0xa1a5e72e'29dbf808'd0fe7e05'869eb825_u128},
615 {Sign::POS, -135, 0xa3884a68'a750cb10'e80a28f4'e1e500d2_u128},
616 {Sign::POS, -135, 0xa58ade36'aeef9f0b'53106415'1ca6e30b_u128},
617 {Sign::POS, -135, 0xa78d7a19'82c4b08f'27c01ffa'8e2e3c4b_u128},
618 {Sign::POS, -135, 0xa9901e11'63cbbbf5'7ba9408d'c857d568_u128},
619 {Sign::POS, -135, 0xab92ca1e'93038d76'104d1e33'31d3b4fa_u128},
620 {Sign::POS, -135, 0xad957e41'516e0158'9343c846'fcdf9137_u128},
621 {Sign::POS, -135, 0xaf780e79'b2514889'3977e89a'ec59bfa2_u128},
622 {Sign::POS, -135, 0xb17ad246'ef3713bc'913d4e3d'c55c3e6e_u128},
623 {Sign::POS, -135, 0xb37d9e2a'7a56b09d'777b52a9'e70d8bcc_u128},
624 {Sign::POS, -135, 0xb5807224'94be0c91'55de916f'd30591de_u128},
625 {Sign::POS, -135, 0xb7834e35'7f7e2600'e79cfb37'be2861e4_u128},
626 {Sign::POS, -135, 0xb986325d'7bab0c89'90983104'd3805389_u128},
627 {Sign::POS, -135, 0xbb68ef9c'254aa378'59e3b2ec'71ce64f4_u128},
628 {Sign::POS, -135, 0xbd6be371'8c77636f'e83183bf'3dd612ef_u128},
629 {Sign::POS, -135, 0xbf6edf5e'c44d9d35'c4e3b0ac'2fd52b7f_u128},
630 };
631
632 // Logarithm range reduction - Step 3:
633 // s(k) = 2^-21 round( 2^21 / (1 + k*2^-21) ) - 1 for k = -69 .. 69
634 // Output range:
635 // [-0x1.012bb800000800114p-22, 0x1p-22 ]
636 constexpr double S3[139] = {
637 0x1.14p-15, 0x1.1p-15, 0x1.0cp-15, 0x1.08p-15, 0x1.04p-15, 0x1p-15,
638 0x1.f8p-16, 0x1.fp-16, 0x1.e8p-16, 0x1.ep-16, 0x1.d8p-16, 0x1.dp-16,
639 0x1.c8p-16, 0x1.cp-16, 0x1.b8p-16, 0x1.bp-16, 0x1.a8p-16, 0x1.ap-16,
640 0x1.98p-16, 0x1.9p-16, 0x1.88p-16, 0x1.8p-16, 0x1.78p-16, 0x1.7p-16,
641 0x1.68p-16, 0x1.6p-16, 0x1.58p-16, 0x1.5p-16, 0x1.48p-16, 0x1.4p-16,
642 0x1.38p-16, 0x1.3p-16, 0x1.28p-16, 0x1.2p-16, 0x1.18p-16, 0x1.1p-16,
643 0x1.08p-16, 0x1p-16, 0x1.fp-17, 0x1.ep-17, 0x1.dp-17, 0x1.cp-17,
644 0x1.bp-17, 0x1.ap-17, 0x1.9p-17, 0x1.8p-17, 0x1.7p-17, 0x1.6p-17,
645 0x1.5p-17, 0x1.4p-17, 0x1.3p-17, 0x1.2p-17, 0x1.1p-17, 0x1p-17,
646 0x1.ep-18, 0x1.cp-18, 0x1.ap-18, 0x1.8p-18, 0x1.6p-18, 0x1.4p-18,
647 0x1.2p-18, 0x1p-18, 0x1.cp-19, 0x1.8p-19, 0x1.4p-19, 0x1p-19,
648 0x1.8p-20, 0x1p-20, 0x1p-21, 0.0, -0x1p-21, -0x1p-20,
649 -0x1.8p-20, -0x1p-19, -0x1.4p-19, -0x1.8p-19, -0x1.cp-19, -0x1p-18,
650 -0x1.2p-18, -0x1.4p-18, -0x1.6p-18, -0x1.8p-18, -0x1.ap-18, -0x1.cp-18,
651 -0x1.ep-18, -0x1p-17, -0x1.1p-17, -0x1.2p-17, -0x1.3p-17, -0x1.4p-17,
652 -0x1.5p-17, -0x1.6p-17, -0x1.7p-17, -0x1.8p-17, -0x1.9p-17, -0x1.ap-17,
653 -0x1.bp-17, -0x1.cp-17, -0x1.dp-17, -0x1.ep-17, -0x1.fp-17, -0x1p-16,
654 -0x1.08p-16, -0x1.1p-16, -0x1.18p-16, -0x1.2p-16, -0x1.28p-16, -0x1.3p-16,
655 -0x1.38p-16, -0x1.4p-16, -0x1.48p-16, -0x1.5p-16, -0x1.58p-16, -0x1.6p-16,
656 -0x1.68p-16, -0x1.7p-16, -0x1.78p-16, -0x1.8p-16, -0x1.88p-16, -0x1.9p-16,
657 -0x1.98p-16, -0x1.ap-16, -0x1.a8p-16, -0x1.bp-16, -0x1.b8p-16, -0x1.cp-16,
658 -0x1.c8p-16, -0x1.dp-16, -0x1.d8p-16, -0x1.ep-16, -0x1.e8p-16, -0x1.fp-16,
659 -0x1.f8p-16, -0x1p-15, -0x1.04p-15, -0x1.08p-15, -0x1.0cp-15, -0x1.1p-15,
660 -0x1.14p-15,
661 };
662
663 // -log(r) for the third step, generated by SageMath with:
664 //
665 // for i in range(-69, 70):
666 // r = 2^-21 * round( 2^21 / (1 + i*2^(-21)) );
667 // s, m, e = RealField(128)(r).log().sign_mantissa_exponent();
668 // print("{Sign::POS," if (s == -1) else "{Sign::NEG,", e, ",
669 // format_hex(m), "},");
670 constexpr Float128 LOG_R3[139] = {
671 {Sign::NEG, -142, 0x89ff6b38'd5de2622'e39d3faf'42340ed7_u128},
672 {Sign::NEG, -142, 0x87ff6f80'ccb40f16'7ff33266'82c02485_u128},
673 {Sign::NEG, -142, 0x85ff73b8'c3cdf731'5caf4fbe'343cf928_u128},
674 {Sign::NEG, -142, 0x83ff77e0'bb2ade79'cdb6e554'348f7fe8_u128},
675 {Sign::NEG, -142, 0x81ff7bf8'b2c9c4f6'0ef009c2'457de25d_u128},
676 {Sign::NEG, -143, 0xffff0001'55535558'8883333c'57b57c74_u128},
677 {Sign::NEG, -143, 0xfbff07f1'45931f44'f32668f3'9c70d183_u128},
678 {Sign::NEG, -143, 0xf7ff0fc1'3650e7bd'459a73c6'a6486fe3_u128},
679 {Sign::NEG, -143, 0xf3ff1771'278aaecd'37b18cca'7dd3a29f_u128},
680 {Sign::NEG, -143, 0xefff1f01'193e7480'513f610d'21bcfc78_u128},
681 {Sign::NEG, -143, 0xebff2671'0b6a38e1'ea190b95'c0690b7b_u128},
682 {Sign::NEG, -143, 0xe7ff2dc0'fe0bfbfd'2a150f64'f0ad1743_u128},
683 {Sign::NEG, -143, 0xe3ff34f0'f121bddd'090b5174'e995e9d1_u128},
684 {Sign::NEG, -143, 0xdfff3c00'e4a97e8c'4ed512b9'b93ea2bf_u128},
685 {Sign::NEG, -143, 0xdbff42f0'd8a13e15'934cea21'7ab794a2_u128},
686 {Sign::NEG, -143, 0xd7ff49c0'cd06fc83'3e4ebe94'8afd2c76_u128},
687 {Sign::NEG, -143, 0xd3ff5070'c1d8b9df'87b7c0f5'bcfee2e1_u128},
688 {Sign::NEG, -143, 0xcfff5700'b7147634'77666622'8cb6371b_u128},
689 {Sign::NEG, -143, 0xcbff5d70'acb8318b'e53a60f3'514db358_u128},
690 {Sign::NEG, -143, 0xc7ff63c0'a2c1ebef'79149c3b'6e57fa86_u128},
691 {Sign::NEG, -143, 0xc3ff69f0'992fa568'aad734c9'8416df2a_u128},
692 {Sign::NEG, -143, 0xbfff7000'8fff5e00'c2657367'9ed28334_u128},
693 {Sign::NEG, -143, 0xbbff75f0'872f15c0'd7a3c6db'6540809f_u128},
694 {Sign::NEG, -143, 0xb7ff7bc0'7ebcccb1'd277bde6'45fb1aad_u128},
695 {Sign::NEG, -143, 0xb3ff8170'76a682dc'6ac80145'a4087793_u128},
696 {Sign::NEG, -143, 0xafff8700'6eea3849'287c4db3'0271e265_u128},
697 {Sign::NEG, -143, 0xabff8c70'6785ed00'637d6de4'2eeb151e_u128},
698 {Sign::NEG, -143, 0xa7ff91c0'6077a10a'43b5348b'6b898a8c_u128},
699 {Sign::NEG, -143, 0xa3ff96f0'59bd546e'c10e7657'978bd7f6_u128},
700 {Sign::NEG, -143, 0x9fff9c00'53550735'a37503f4'57310e59_u128},
701 {Sign::NEG, -143, 0x9bffa0f0'4d3cb966'82d5a40a'3aa022ff_u128},
702 {Sign::NEG, -143, 0x97ffa5c0'47726b08'c71e0d3e'e3df5f4d_u128},
703 {Sign::NEG, -143, 0x93ffaa70'41f41c23'a83ce035'2bdbd79b_u128},
704 {Sign::NEG, -143, 0x8fffaf00'3cbfccbe'2e21a18d'4680e8e4_u128},
705 {Sign::NEG, -143, 0x8bffb370'37d37cdf'30bcb3e4'e5dfbd28_u128},
706 {Sign::NEG, -143, 0x87ffb7c0'332d2c8d'57ff51d7'5c66d64a_u128},
707 {Sign::NEG, -143, 0x83ffbbf0'2ecadbcf'1bdb87fd'be299f43_u128},
708 {Sign::NEG, -144, 0xffff8000'55551555'88885dde'02700703_u128},
709 {Sign::NEG, -144, 0xf7ff87e0'4d94724c'd259ca80'3a0c1870_u128},
710 {Sign::NEG, -144, 0xefff8f80'464fce8f'e5141308'51c7070a_u128},
711 {Sign::NEG, -144, 0xe7ff96e0'3f832a2a'30a16898'f3073a64_u128},
712 {Sign::NEG, -144, 0xdfff9e00'392a8526'c4ed6451'7b2949ce_u128},
713 {Sign::NEG, -144, 0xd7ffa4e0'3341df90'51e4fb4e'32cf6350_u128},
714 {Sign::NEG, -144, 0xcfffab80'2dc53971'277672a8'8350bcce_u128},
715 {Sign::NEG, -144, 0xc7ffb1e0'28b092d3'35915377'2a490f06_u128},
716 {Sign::NEG, -144, 0xbfffb800'23ffebc0'0c265ece'6b481a0e_u128},
717 {Sign::NEG, -144, 0xb7ffbde0'1faf4440'db2781c0'3fa132f6_u128},
718 {Sign::NEG, -144, 0xafffc380'1bba9c5e'7287c95c'845ada33_u128},
719 {Sign::NEG, -144, 0xa7ffc8e0'181df421'423b56b1'263e5a77_u128},
720 {Sign::NEG, -144, 0x9fffce00'14d54b91'5a3752ca'4c076fa3_u128},
721 {Sign::NEG, -144, 0x97ffd2e0'11dca2b6'6a71e2b2'7eb3f573_u128},
722 {Sign::NEG, -144, 0x8fffd780'0f2ff997'c2e21b72'cff39d8f_u128},
723 {Sign::NEG, -144, 0x87ffdbe0'0ccb503c'537ff612'feb7ac9e_u128},
724 {Sign::NEG, -145, 0xffffc000'15554d55'58888733'33c57c18_u128},
725 {Sign::NEG, -145, 0xefffc7c0'1193f9d1'fa514218'42311c42_u128},
726 {Sign::NEG, -145, 0xdfffcf00'0e4aa5fa'2c4ed6de'475b942c_u128},
727 {Sign::NEG, -145, 0xcfffd5c0'0b7151d8'ce77678c'bb6fcb88_u128},
728 {Sign::NEG, -145, 0xbfffdc00'08fffd78'00c26629'a679ed3b_u128},
729 {Sign::NEG, -145, 0xafffe1c0'06eea8e1'23287cb9'd3072728_u128},
730 {Sign::NEG, -145, 0x9fffe700'0535541c'd5a37540'fd057315_u128},
731 {Sign::NEG, -145, 0x8fffebc0'03cbff32'f82e21c1'fce36810_u128},
732 {Sign::NEG, -146, 0xffffe000'05555455'5588887d'dde02702_u128},
733 {Sign::NEG, -146, 0xdfffe780'0392aa14'9ac4ed72'adf5b295_u128},
734 {Sign::NEG, -146, 0xbfffee00'023fffaf'000c2664'8066b482_u128},
735 {Sign::NEG, -146, 0x9ffff380'014d552e'455a3754'b292c077_u128},
736 {Sign::NEG, -147, 0xfffff000'01555535'55588888'33333c58_u128},
737 {Sign::NEG, -147, 0xbffff700'008ffff5'e000c266'5736679f_u128},
738 {Sign::NEG, -148, 0xfffff800'00555551'55558888'85ddde02_u128},
739 {Sign::NEG, -149, 0xfffffc00'00155554'd5555888'88733334_u128},
740 {Sign::POS, 0, 0_u128},
741 {Sign::POS, -148, 0x80000200'000aaaaa'eaaaac44'444eeeef_u128},
742 {Sign::POS, -147, 0x80000400'002aaaac'aaaac444'459999ac_u128},
743 {Sign::POS, -147, 0xc0000900'0090000a'2000c266'7596679f_u128},
744 {Sign::POS, -146, 0x80000800'00aaaaba'aaac4444'6eeef381_u128},
745 {Sign::POS, -146, 0xa0000c80'014d557c'655a3755'f81815cc_u128},
746 {Sign::POS, -146, 0xc0001200'02400051'000c2668'4c66b482_u128},
747 {Sign::POS, -146, 0xe0001880'0392ab40'bac4ed7c'40fb07eb_u128},
748 {Sign::POS, -145, 0x80001000'02aaab2a'aac44449'999abe2c_u128},
749 {Sign::POS, -145, 0x90001440'03cc00cd'082e21d7'9cbb6812_u128},
750 {Sign::POS, -145, 0xa0001900'0535568d'd5a37569'adb01dc3_u128},
751 {Sign::POS, -145, 0xb0001e40'06eeac74'33287d01'e8c9d1d9_u128},
752 {Sign::POS, -145, 0xc0002400'09000288'00c266a3'2679ed48_u128},
753 {Sign::POS, -145, 0xd0002a40'0b7158d1'de776851'22b2764b_u128},
754 {Sign::POS, -145, 0xe0003100'0e4aaf5b'2c4ed810'a8063f03_u128},
755 {Sign::POS, -145, 0xf0003840'1194062e'0a5143e7'be891c8f_u128},
756 {Sign::POS, -144, 0x80002000'0aaaaeaa'ac4444ee'ef3813a1_u128},
757 {Sign::POS, -144, 0x88002420'0ccb5a6e'5b7ff7fe'1339025b_u128},
758 {Sign::POS, -144, 0x90002880'0f300668'42e21e26'caf39e33_u128},
759 {Sign::POS, -144, 0x98002d20'11dcb29e'f271e66f'a5554bc6_u128},
760 {Sign::POS, -144, 0xa0003200'14d55f19'5a3757e0'615cc676_u128},
761 {Sign::POS, -144, 0xa8003720'181e0bde'ca3b5d82'10ca5cab_u128},
762 {Sign::POS, -144, 0xb0003c80'1bbab8f6'f287d25f'3cb032bb_u128},
763 {Sign::POS, -144, 0xb8004220'1faf6669'e3278d84'0be28cdb_u128},
764 {Sign::POS, -144, 0xc0004800'24001440'0c266dfe'6b482076_u128},
765 {Sign::POS, -144, 0xc8004e20'28b0c282'3d9166de'380a6d3d_u128},
766 {Sign::POS, -144, 0xd0005480'2dc57139'a7768b35'6ba61e4b_u128},
767 {Sign::POS, -144, 0xd8005b20'3342206f'd9e51a18'49db73c1_u128},
768 {Sign::POS, -144, 0xe0006200'392ad02e'c4ed8a9d'907eb521_u128},
769 {Sign::POS, -144, 0xe8006920'3f838080'b8a197de'a928acd7_u128},
770 {Sign::POS, -144, 0xf0007080'46503170'65144cf7'dcc72d3b_u128},
771 {Sign::POS, -144, 0xf8007820'4d94e308'da5a1108'890d9f6a_u128},
772 {Sign::POS, -143, 0x80004000'2aaacaaa'c4445999'abe2ce2c_u128},
773 {Sign::POS, -143, 0x84004410'2ecb2431'1fdbbb4f'3bffc832_u128},
774 {Sign::POS, -143, 0x88004840'332d7e1d'97ff8f39'ec91b4ee_u128},
775 {Sign::POS, -143, 0x8c004c90'37d3d876'74bcfcf0'b3f0a95d_u128},
776 {Sign::POS, -143, 0x90005100'3cc03342'2e21f80c'a6813aff_u128},
777 {Sign::POS, -143, 0x94005590'41f48e87'6c3d4629'170ce87f_u128},
778 {Sign::POS, -143, 0x98005a40'4772ea4d'071e84e3'b80a8881_u128},
779 {Sign::POS, -143, 0x9c005f10'4d3d469a'06d62fdc'bdd6bec3_u128},
780 {Sign::POS, -143, 0xa0006400'5355a375'a375a6b7'01dc77c0_u128},
781 {Sign::POS, -143, 0xa4006910'59be00e7'450f3318'26ad6b05_u128},
782 {Sign::POS, -143, 0xa8006e40'60785ef6'83b60ea8'bd0aa459_u128},
783 {Sign::POS, -143, 0xac007390'6786bdab'277e6914'69dd13f5_u128},
784 {Sign::POS, -143, 0xb0007900'6eeb1d0d'287d6e0a'0d1e25eb_u128},
785 {Sign::POS, -143, 0xb4007e90'76a77d24'aec94b3b'e9b060f5_u128},
786 {Sign::POS, -143, 0xb8008440'7ebdddfa'1279365f'ce280cce_u128},
787 {Sign::POS, -143, 0xbc008a10'87303f95'dba5732f'3e83e04a_u128},
788 {Sign::POS, -143, 0xc0009000'9000a200'c2675967'9ed5b754_u128},
789 {Sign::POS, -143, 0xc4009610'99310543'aed95aca'5edb5109_u128},
790 {Sign::POS, -143, 0xc8009c40'a2c36967'b917091d'2687160f_u128},
791 {Sign::POS, -143, 0xcc00a290'acb9ce76'293d1c2a'0378e75d_u128},
792 {Sign::POS, -143, 0xd000a900'b7163478'776977bf'9766f5a7_u128},
793 {Sign::POS, -143, 0xd400af90'c1da9b78'4bbb31b1'4776a18b_u128},
794 {Sign::POS, -143, 0xd800b640'cd09037f'7e5297d7'6c8564ba_u128},
795 {Sign::POS, -143, 0xdc00bd10'd8a36c98'1751360f'8461c447_u128},
796 {Sign::POS, -143, 0xe000c400'e4abd6cc'4ed9dc3c'63f44c41_u128},
797 {Sign::POS, -143, 0xe400cb10'f1244226'8d10a446'6a5894d5_u128},
798 {Sign::POS, -143, 0xe800d240'fe0eaeb1'6a1af81b'b4e6510e_u128},
799 {Sign::POS, -143, 0xec00d991'0b6d1c77'ae1f97b0'542a677a_u128},
800 {Sign::POS, -143, 0xf000e101'19418b84'51469efe'81d014cc_u128},
801 {Sign::POS, -143, 0xf400e891'278dfbe2'7bb98c06'd77a18b4_u128},
802 {Sign::POS, -143, 0xf800f041'36546d9d'85a344d0'868bed17_u128},
803 {Sign::POS, -143, 0xfc00f811'4596e0c0'f7301d69'90e307cc_u128},
804 {Sign::POS, -142, 0x80008000'aaabaaac'4446eef3'8140138f_u128},
805 {Sign::POS, -142, 0x82008408'b2cbe5b8'10f5e432'96105497_u128},
806 {Sign::POS, -142, 0x84008820'bb2d2189'edbd4f83'ef63f730_u128},
807 {Sign::POS, -142, 0x86008c48'c3d05e27'feb654fd'541c638e_u128},
808 {Sign::POS, -142, 0x88009080'ccb69b98'7ffadeb8'882f7674_u128},
809 {Sign::POS, -142, 0x8a0094c8'd5e0d9e1'c5a59fd3'6bd44397_u128},
810 };
811
812 // Minimax polynomial generated by Sollya with:
813 // > P = fpminimax((log(1 + x) - x)/x^2, 3, [|1, 128...|],
814 // [-0x1.01928p-22 , 0x1p-22]);
815 // > P;
816 // > dirtyinfnorm(log(1 + x)/x - 1 - x*P, [-0x1.01928p-22 , 0x1p-22]);
817 // 0x1.ce1e...p-116
818 constexpr Float128 BIG_COEFFS[4]{
819 {Sign::POS, -130, 0xccccccd7'4818e397'7ed78465'd460315b_u128},
820 {Sign::NEG, -129, 0x80000000'000478b0'c6388a23'871ce156_u128},
821 {Sign::POS, -129, 0xaaaaaaaa'aaaaaaaa'aa807bd8'67763262_u128},
822 {Sign::NEG, -128, 0x80000000'00000000'00000000'00000000_u128},
823 };
824
log1p_accurate(int e_x,int index,fputil::DoubleDouble m_x)825 [[maybe_unused]] LIBC_INLINE double log1p_accurate(int e_x, int index,
826 fputil::DoubleDouble m_x) {
827 Float128 e_x_f128(static_cast<float>(e_x));
828 Float128 sum = fputil::quick_mul(LOG_2, e_x_f128);
829 sum = fputil::quick_add(sum, LOG_R1[index]);
830
831 // fputil::DoubleDouble v4;
832 Float128 v = fputil::quick_add(Float128(m_x.hi), Float128(m_x.lo));
833
834 // Skip 2nd range reduction step if |m_x| <= 2^-15.
835 if (m_x.hi > 0x1p-15 || m_x.hi < -0x1p-15) {
836 // Range reduction - Step 2.
837 // For k such that: k * 2^-14 - 2^-15 <= m_x.hi < k * 2^-14 + 2^-15,
838 // Let s_k = 2^-18 * round( 2^18 / (1 + k*2^-14) ) - 1
839 // Then the 2nd reduced argument is:
840 // (1 + s_k) * (1 + m_x) - 1 =
841 // = s_k + m_x + s_k * m_x
842 // Output range:
843 // -0x1.1037c00000040271p-15 <= v2.hi + v2.lo <= 0x1.108480000008096cp-15
844 int idx2 = static_cast<int>(0x1p14 * (m_x.hi + (91 * 0x1p-14 + 0x1p-15)));
845 sum = fputil::quick_add(sum, LOG_R2[idx2]);
846 Float128 s2 = Float128(S2[idx2]);
847 v = fputil::quick_add(fputil::quick_add(v, s2), fputil::quick_mul(v, s2));
848 }
849
850 // Skip 3rd range reduction step if |v| <= 2^-22.
851 if (v.exponent > -150) {
852 // Range reduction - Step 3.
853 // For k such that: k * 2^-21 - 2^-22 <= v2.hi < k * 2^-21 + 2^-22,
854 // Let s_k = 2^-21 * round( 2^21 / (1 + k*2^-21) ) - 1
855 // Then the 3rd reduced argument is:
856 // v3.hi + v3.lo ~ (1 + s_k) * (1 + v2.hi + v2.lo) - 1
857 // Output range:
858 // -0x1.012bb800000800114p-22 <= v3.hi + v3.lo <= 0x1p-22
859 int idx3 =
860 static_cast<int>(0x1p21 * (double(v) + (69 * 0x1p-21 + 0x1p-22)));
861 sum = fputil::quick_add(sum, LOG_R3[idx3]);
862 Float128 s3 = Float128(S3[idx3]);
863 v = fputil::quick_add(fputil::quick_add(v, s3), fputil::quick_mul(v, s3));
864 }
865
866 // Polynomial approximation
867 Float128 p = fputil::quick_mul(v, BIG_COEFFS[0]);
868 p = fputil::quick_mul(v, fputil::quick_add(p, BIG_COEFFS[1]));
869 p = fputil::quick_mul(v, fputil::quick_add(p, BIG_COEFFS[2]));
870 p = fputil::quick_mul(v, fputil::quick_add(p, BIG_COEFFS[3]));
871 p = fputil::quick_add(v, fputil::quick_mul(v, p));
872
873 Float128 r = fputil::quick_add(sum, p);
874
875 return static_cast<double>(r);
876 }
877
878 } // namespace
879
880 LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
881 using FPBits_t = typename fputil::FPBits<double>;
882
883 constexpr int EXP_BIAS = FPBits_t::EXP_BIAS;
884 constexpr int FRACTION_LEN = FPBits_t::FRACTION_LEN;
885 FPBits_t xbits(x);
886 uint64_t x_u = xbits.uintval();
887
888 fputil::DoubleDouble x_dd{0.0, 0.0};
889
890 uint16_t x_exp = xbits.get_biased_exponent();
891
892 if (x_exp >= EXP_BIAS) {
893 // |x| >= 1
894 if (LIBC_UNLIKELY(x_u >= 0x4650'0000'0000'0000ULL)) {
895 // x >= 2^102 or x is negative, inf, or NaN
896 if (LIBC_UNLIKELY(x_u > FPBits_t::max_normal().uintval())) {
897 // x <= -1.0 or x is Inf or NaN
898 if (x_u == 0xbff0'0000'0000'0000ULL) {
899 // x = -1.0
900 fputil::set_errno_if_required(ERANGE);
901 fputil::raise_except_if_required(FE_DIVBYZERO);
902 return FPBits_t::inf(Sign::NEG).get_val();
903 }
904 if (xbits.is_neg() && !xbits.is_nan()) {
905 // x < -1.0
906 fputil::set_errno_if_required(EDOM);
907 fputil::raise_except_if_required(FE_INVALID);
908 return FPBits_t::quiet_nan().get_val();
909 }
910 // x is +Inf or NaN
911 return x;
912 }
913 x_dd.hi = x;
914 } else {
915 x_dd = fputil::exact_add(x, 1.0);
916 }
917 } else {
918 // |x| < 1
919 if (LIBC_UNLIKELY(xbits.get_biased_exponent() <
920 EXP_BIAS - FRACTION_LEN - 1)) {
921 // Quick return when |x| < 2^-53.
922 // Since log(1 + x) = x - x^2/2 + x^3/3 - ...,
923 // for |x| < 2^-53,
924 // x > log(1 + x) > x - x^2 > x(1 - 2^-54) > x - ulp(x)/2
925 // Thus,
926 // log(1 + x) = nextafter(x, -inf) for FE_DOWNWARD, or
927 // FE_TOWARDZERO and x > 0,
928 // = x otherwise.
929 if (x == 0.0)
930 return x + x; // Handle FTZ/DAZ correctly.
931
932 volatile float tp = 1.0f;
933 volatile float tn = -1.0f;
934 bool rdp = (tp - 0x1p-25f != tp);
935 bool rdn = (tn - 0x1p-24f != tn);
936
937 if (x > 0 && rdp) {
938 return FPBits_t(x_u - 1).get_val();
939 }
940
941 if (x < 0 && rdn) {
942 return FPBits_t(x_u + 1).get_val();
943 }
944
945 return (x + x == 0.0) ? x + x : x;
946 }
947 x_dd = fputil::exact_add(1.0, x);
948 }
949
950 // At this point, x_dd is the exact sum of 1 + x:
951 // x_dd.hi + x_dd.lo = x + 1.0 exactly.
952 // |x_dd.hi| >= 2^-54
953 // |x_dd.lo| < ulp(x_dd.hi)
954
955 FPBits_t xhi_bits(x_dd.hi);
956 uint64_t xhi_frac = xhi_bits.get_mantissa();
957 x_u = xhi_bits.uintval();
958 // Range reduction:
959 // Find k such that |x_hi - k * 2^-7| <= 2^-8.
960 int idx = static_cast<int>((xhi_frac + (1ULL << (FRACTION_LEN - 8))) >>
961 (FRACTION_LEN - 7));
962 int x_e = xhi_bits.get_exponent() + (idx >> 7);
963 double e_x = static_cast<double>(x_e);
964
965 // hi is exact
966 // ulp(hi) = ulp(LOG_2_HI) = ulp(LOG_R1_DD[idx].hi) = 2^-43
967 double hi = fputil::multiply_add(e_x, LOG_2_HI, LOG_R1_DD[idx].hi);
968 // lo errors < |e_x| * ulp(LOG_2_LO) + ulp(LOG_R1[idx].lo)
969 // <= 2^11 * 2^(-43-53) = 2^-85
970 double lo = fputil::multiply_add(e_x, LOG_2_LO, LOG_R1_DD[idx].lo);
971
972 // Error bound of e_x * log(2) - log(r1)
973 constexpr double ERR_HI[2] = {0x1.0p-85, 0.0};
974 double err_hi = ERR_HI[hi == 0.0];
975
976 // Scale x_dd by 2^(-xh_bits.get_exponent()).
977 int64_t s_u = static_cast<int64_t>(x_u & FPBits_t::EXP_MASK) -
978 (static_cast<int64_t>(EXP_BIAS) << FRACTION_LEN);
979 // Normalize arguments:
980 // 1 <= m_dd.hi < 2
981 // |m_dd.lo| < 2^-52.
982 // This is exact.
983 uint64_t m_hi = FPBits_t::one().uintval() | xhi_frac;
984
985 uint64_t m_lo =
986 FPBits_t(x_dd.lo).abs().get_val() > x_dd.hi * 0x1.0p-127
987 ? static_cast<uint64_t>(cpp::bit_cast<int64_t>(x_dd.lo) - s_u)
988 : 0;
989
990 fputil::DoubleDouble m_dd{FPBits_t(m_lo).get_val(), FPBits_t(m_hi).get_val()};
991
992 // Perform range reduction:
993 // r * m - 1 = r * (m_dd.hi + m_dd.lo) - 1
994 // = (r * m_dd.hi - 1) + r * m_dd.lo
995 // = v_hi + (v_lo.hi + v_lo.lo)
996 // where:
997 // v_hi = r * m_dd.hi - 1 (exact)
998 // v_lo.hi + v_lo.lo = r * m_dd.lo (exact)
999 // Bounds on the values:
1000 // -0x1.69000000000edp-8 < r * m - 1 < 0x1.7f00000000081p-8
1001 // |v_lo.hi| <= |r| * |m_dd.lo| < 2^-52
1002 // |v_lo.lo| < ulp(v_lo.hi) <= 2^(-52 - 53) = 2^(-105)
1003 double r = R1[idx];
1004 double v_hi;
1005 fputil::DoubleDouble v_lo = fputil::exact_mult(m_dd.lo, r);
1006
1007 // Perform exact range reduction
1008 #ifdef LIBC_TARGET_CPU_HAS_FMA
1009 v_hi = fputil::multiply_add(r, m_dd.hi, -1.0); // Exact.
1010 #else
1011 // c = 1 + idx * 2^-7.
1012 double c = FPBits_t((static_cast<uint64_t>(idx) << (FRACTION_LEN - 7)) +
1013 uint64_t(0x3FF0'0000'0000'0000ULL))
1014 .get_val();
1015 v_hi = fputil::multiply_add(r, m_dd.hi - c, RCM1[idx]); // Exact
1016 #endif // LIBC_TARGET_CPU_HAS_FMA
1017
1018 // Range reduction output:
1019 // -0x1.69000000000edp-8 < v_hi + v_lo < 0x1.7f00000000081p-8
1020 // |v_dd.lo| < ulp(v_dd.hi) <= 2^(-7 - 53) = 2^-60
1021 fputil::DoubleDouble v_dd = fputil::exact_add(v_hi, v_lo.hi);
1022 v_dd.lo += v_lo.lo;
1023
1024 // Exact sum:
1025 // r1.hi + r1.lo = e_x * log(2)_hi - log(r)_hi + u
1026 fputil::DoubleDouble r1 = fputil::exact_add(hi, v_dd.hi);
1027
1028 // Overall error is bounded by:
1029 // C * ulp(v_sq) + err_hi
1030 double v_sq = v_dd.hi * v_dd.hi;
1031 double p0 = fputil::multiply_add(v_dd.hi, P_COEFFS[1], P_COEFFS[0]);
1032 double p1 = fputil::multiply_add(v_dd.hi, P_COEFFS[3], P_COEFFS[2]);
1033 double p2 = fputil::multiply_add(v_dd.hi, P_COEFFS[5], P_COEFFS[4]);
1034 double p = fputil::polyeval(v_sq, (v_dd.lo + r1.lo) + lo, p0, p1, p2);
1035
1036 double err = fputil::multiply_add(v_sq, P_ERR, err_hi);
1037
1038 double left = r1.hi + (p - err);
1039 double right = r1.hi + (p + err);
1040
1041 // Ziv's test to see if fast pass is accurate enough.
1042 if (left == right)
1043 return left;
1044
1045 return log1p_accurate(x_e, idx, v_dd);
1046 }
1047
1048 } // namespace LIBC_NAMESPACE_DECL
1049