[libc-commits] [libc] b91e78d - [libc][math] Implement double precision log1p correctly rounded to all rounding modes.
Tue Ly via libc-commits
libc-commits at lists.llvm.org
Tue May 23 08:04:21 PDT 2023
Author: Tue Ly
Date: 2023-05-23T11:04:04-04:00
New Revision: b91e78da37800c03542141135a6628f320e974ed
URL: https://github.com/llvm/llvm-project/commit/b91e78da37800c03542141135a6628f320e974ed
DIFF: https://github.com/llvm/llvm-project/commit/b91e78da37800c03542141135a6628f320e974ed.diff
LOG: [libc][math] Implement double precision log1p correctly rounded to all rounding modes.
Implement double precision log1p function correctly rounded to all
rounding modes.
**Performance**
- For `0.5 <= x <= 2`, the fast pass hitting rate is about 99.93%.
- Benchmarks with `./perf.sh` tool from the CORE-MATH project, unit is (CPU clocks / call).
- Reciprocal throughput from CORE-MATH's perf tool on Ryzen 5900X:
```
$ ./perf.sh log1p
GNU libc version: 2.35
GNU libc release: stable
-- CORE-MATH reciprocal throughput -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 39.792 + 1.011 clc/call; Median-Min = 0.940 clc/call; Max = 41.373 clc/call;
-- CORE-MATH reciprocal throughput -- without FMA (-march=x86-64-v2)
[####################] 100 %
Ntrial = 20 ; Min = 87.285 + 1.135 clc/call; Median-Min = 1.299 clc/call; Max = 89.715 clc/call;
-- System LIBC reciprocal throughput --
[####################] 100 %
Ntrial = 20 ; Min = 20.666 + 0.123 clc/call; Median-Min = 0.125 clc/call; Max = 20.828 clc/call;
-- LIBC reciprocal throughput -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 20.928 + 0.771 clc/call; Median-Min = 0.725 clc/call; Max = 22.767 clc/call;
-- LIBC reciprocal throughput -- without FMA
[####################] 100 %
Ntrial = 20 ; Min = 31.461 + 0.528 clc/call; Median-Min = 0.602 clc/call; Max = 36.809 clc/call;
```
- Latency from CORE-MATH's perf tool on Ryzen 5900X:
```
$ ./perf.sh log1p --latency
GNU libc version: 2.35
GNU libc release: stable
-- CORE-MATH latency -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 77.875 + 0.062 clc/call; Median-Min = 0.051 clc/call; Max = 78.003 clc/call;
-- CORE-MATH latency -- without FMA (-march=x86-64-v2)
[####################] 100 %
Ntrial = 20 ; Min = 101.958 + 1.202 clc/call; Median-Min = 1.325 clc/call; Max = 104.452 clc/call;
-- System LIBC latency --
[####################] 100 %
Ntrial = 20 ; Min = 60.581 + 1.443 clc/call; Median-Min = 1.611 clc/call; Max = 62.285 clc/call;
-- LIBC latency -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 48.817 + 1.108 clc/call; Median-Min = 1.300 clc/call; Max = 50.282 clc/call;
-- LIBC latency -- without FMA
[####################] 100 %
Ntrial = 20 ; Min = 61.121 + 0.599 clc/call; Median-Min = 0.761 clc/call; Max = 62.020 clc/call;
```
- Accurate pass latency:
```
$ ./perf.sh log1p --latency --simple_stat
GNU libc version: 2.35
GNU libc release: stable
-- CORE-MATH latency -- with FMA
760.444
-- CORE-MATH latency -- without FMA (-march=x86-64-v2)
827.880
-- LIBC latency -- with FMA
711.837
-- LIBC latency -- without FMA
764.317
```
Reviewed By: zimmermann6
Differential Revision: https://reviews.llvm.org/D151049
Added:
libc/src/math/generic/log1p.cpp
libc/src/math/log1p.h
libc/test/src/math/log1p_test.cpp
Modified:
libc/config/darwin/arm/entrypoints.txt
libc/config/linux/aarch64/entrypoints.txt
libc/config/linux/x86_64/entrypoints.txt
libc/config/windows/entrypoints.txt
libc/spec/stdc.td
libc/src/math/CMakeLists.txt
libc/src/math/generic/CMakeLists.txt
libc/test/src/math/CMakeLists.txt
Removed:
################################################################################
diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt
index d5e80ef3a5ad2..0c139c3a04af8 100644
--- a/libc/config/darwin/arm/entrypoints.txt
+++ b/libc/config/darwin/arm/entrypoints.txt
@@ -176,6 +176,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexpl
libc.src.math.log10
libc.src.math.log10f
+ libc.src.math.log1p
libc.src.math.log1pf
libc.src.math.log2
libc.src.math.log2f
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index c3b2fb919d187..2020a45825563 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -287,6 +287,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexpl
libc.src.math.log10
libc.src.math.log10f
+ libc.src.math.log1p
libc.src.math.log1pf
libc.src.math.log2
libc.src.math.log2f
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 31f8b58deee08..546e792e81e8d 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -292,6 +292,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.llroundl
libc.src.math.log10
libc.src.math.log10f
+ libc.src.math.log1p
libc.src.math.log1pf
libc.src.math.log2
libc.src.math.log2f
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 6df29390fffb5..3a980a0e5045e 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -169,6 +169,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.llroundl
libc.src.math.log10
libc.src.math.log10f
+ libc.src.math.log1p
libc.src.math.log1pf
libc.src.math.log2
libc.src.math.log2f
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index 8b091b87e749a..be052547c2fac 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -408,6 +408,7 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"log10", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"log10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+ FunctionSpec<"log1p", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"log1pf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"log2", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 79d343bf79708..c47584039d5dd 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -114,6 +114,7 @@ add_math_entrypoint_object(ldexpl)
add_math_entrypoint_object(log10)
add_math_entrypoint_object(log10f)
+add_math_entrypoint_object(log1p)
add_math_entrypoint_object(log1pf)
add_math_entrypoint_object(log2)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 3d02851f31e6d..3c2ac5a6f1f34 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -814,6 +814,26 @@ add_entrypoint_object(
-O3
)
+add_entrypoint_object(
+ log1p
+ SRCS
+ log1p.cpp
+ HDRS
+ ../log1p.h
+ DEPENDS
+ .common_constants
+ .log_range_reduction
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.dyadic_float
+ libc.src.__support.macros.optimization
+ COMPILE_OPTIONS
+ -O3
+)
+
add_entrypoint_object(
log1pf
SRCS
diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
new file mode 100644
index 0000000000000..1a944445314e7
--- /dev/null
+++ b/libc/src/math/generic/log1p.cpp
@@ -0,0 +1,1038 @@
+//===-- Double-precision log1p(x) function --------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/log1p.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+#include "common_constants.h"
+#include "log_range_reduction.h"
+
+namespace __llvm_libc {
+
+// 128-bit precision dyadic floating point numbers.
+using Float128 = typename fputil::DyadicFloat<128>;
+using MType = typename Float128::MantissaType;
+
+namespace {
+
+// Extra errors from P is from using x^2 to reduce evaluation latency.
+constexpr double P_ERR = 0x1.0p-50;
+
+// log(2) with 128-bit prepcision generated by SageMath with:
+// sage: (s, m, e) = RealField(128)(2).log().sign_mantissa_exponent();
+// sage: print("MType({", hex(m % 2^64), ",", hex((m >> 64) % 2^64), "})");
+const Float128 LOG_2(/*sign=*/false, /*exponent=*/-128, /*mantissa=*/
+ MType({0xc9e3b39803f2f6af, 0xb17217f7d1cf79ab}));
+
+// R1[i] = 2^-8 * nearestint( 2^8 / (1 + i * 2^-7) )
+constexpr double R1[129] = {
+ 0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.eap-1,
+ 0x1.e6p-1, 0x1.e2p-1, 0x1.dep-1, 0x1.dap-1, 0x1.d8p-1, 0x1.d4p-1, 0x1.dp-1,
+ 0x1.cep-1, 0x1.cap-1, 0x1.c8p-1, 0x1.c4p-1, 0x1.cp-1, 0x1.bep-1, 0x1.bap-1,
+ 0x1.b8p-1, 0x1.b4p-1, 0x1.b2p-1, 0x1.bp-1, 0x1.acp-1, 0x1.aap-1, 0x1.a6p-1,
+ 0x1.a4p-1, 0x1.a2p-1, 0x1.9ep-1, 0x1.9cp-1, 0x1.9ap-1, 0x1.98p-1, 0x1.94p-1,
+ 0x1.92p-1, 0x1.9p-1, 0x1.8ep-1, 0x1.8ap-1, 0x1.88p-1, 0x1.86p-1, 0x1.84p-1,
+ 0x1.82p-1, 0x1.8p-1, 0x1.7ep-1, 0x1.7ap-1, 0x1.78p-1, 0x1.76p-1, 0x1.74p-1,
+ 0x1.72p-1, 0x1.7p-1, 0x1.6ep-1, 0x1.6cp-1, 0x1.6ap-1, 0x1.68p-1, 0x1.66p-1,
+ 0x1.64p-1, 0x1.62p-1, 0x1.6p-1, 0x1.5ep-1, 0x1.5cp-1, 0x1.5ap-1, 0x1.58p-1,
+ 0x1.58p-1, 0x1.56p-1, 0x1.54p-1, 0x1.52p-1, 0x1.5p-1, 0x1.4ep-1, 0x1.4cp-1,
+ 0x1.4ap-1, 0x1.4ap-1, 0x1.48p-1, 0x1.46p-1, 0x1.44p-1, 0x1.42p-1, 0x1.42p-1,
+ 0x1.4p-1, 0x1.3ep-1, 0x1.3cp-1, 0x1.3cp-1, 0x1.3ap-1, 0x1.38p-1, 0x1.36p-1,
+ 0x1.36p-1, 0x1.34p-1, 0x1.32p-1, 0x1.3p-1, 0x1.3p-1, 0x1.2ep-1, 0x1.2cp-1,
+ 0x1.2cp-1, 0x1.2ap-1, 0x1.28p-1, 0x1.28p-1, 0x1.26p-1, 0x1.24p-1, 0x1.24p-1,
+ 0x1.22p-1, 0x1.2p-1, 0x1.2p-1, 0x1.1ep-1, 0x1.1cp-1, 0x1.1cp-1, 0x1.1ap-1,
+ 0x1.1ap-1, 0x1.18p-1, 0x1.16p-1, 0x1.16p-1, 0x1.14p-1, 0x1.14p-1, 0x1.12p-1,
+ 0x1.12p-1, 0x1.1p-1, 0x1.0ep-1, 0x1.0ep-1, 0x1.0cp-1, 0x1.0cp-1, 0x1.0ap-1,
+ 0x1.0ap-1, 0x1.08p-1, 0x1.08p-1, 0x1.06p-1, 0x1.06p-1, 0x1.04p-1, 0x1.04p-1,
+ 0x1.02p-1, 0x1.02p-1, 0x1p-1,
+};
+
+// Extra constants for exact range reduction when FMA instructions are not
+// available:
+// r * c - 1 for r = 2^-8 * nearestint( 2^8 / (1 + i * 2^-7))
+// and c = 1 + i * 2^-7
+// with i = 0..128.
+[[maybe_unused]] constexpr double RCM1[129] = {
+ 0.0, -0x1p-14, -0x1p-12, -0x1.2p-11, -0x1p-10, -0x1.9p-10,
+ 0x1.fp-10, 0x1.28p-10, 0x1p-12, -0x1.9p-11, -0x1.fp-10, 0x1.2p-10,
+ -0x1p-12, -0x1.cp-10, 0x1.1p-10, -0x1.5p-11, 0x1p-9, 0x1p-14,
+ -0x1p-9, 0x1.ap-12, -0x1.ep-10, 0x1.8p-12, -0x1.1p-9, -0x1p-15,
+ 0x1p-9, -0x1.ap-11, 0x1.1p-10, -0x1.f8p-10, -0x1p-12, 0x1.68p-10,
+ -0x1.fp-10, -0x1.cp-12, 0x1p-10, 0x1.3p-9, -0x1.6p-10, -0x1.4p-13,
+ 0x1p-10, 0x1.0cp-9, -0x1.08p-9, -0x1.2p-10, -0x1p-12, 0x1.2p-11,
+ 0x1.5p-10, 0x1p-9, 0x1.5p-9, -0x1.1cp-9, -0x1.cp-10, -0x1.58p-10,
+ -0x1p-10, -0x1.7p-11, -0x1p-11, -0x1.6p-12, -0x1p-12, -0x1.cp-13,
+ -0x1p-12, -0x1.6p-12, -0x1p-11, -0x1.7p-11, -0x1p-10, -0x1.58p-10,
+ -0x1.cp-10, -0x1.1cp-9, -0x1.6p-9, 0x1.5p-9, 0x1p-9, 0x1.5p-10,
+ 0x1.2p-11, -0x1p-12, -0x1.2p-10, -0x1.08p-9, -0x1.88p-9, 0x1.0cp-9,
+ 0x1p-10, -0x1.4p-13, -0x1.6p-10, -0x1.54p-9, 0x1.3p-9, 0x1p-10,
+ -0x1.cp-12, -0x1.fp-10, 0x1.8p-9, 0x1.68p-10, -0x1p-12, -0x1.f8p-10,
+ 0x1.7p-9, 0x1.1p-10, -0x1.ap-11, -0x1.6p-9, 0x1p-9, -0x1p-15,
+ -0x1.1p-9, 0x1.48p-9, 0x1.8p-12, -0x1.ep-10, 0x1.6p-9, 0x1.ap-12,
+ -0x1p-9, 0x1.48p-9, 0x1p-14, -0x1.4p-9, 0x1p-9, -0x1.5p-11,
+ -0x1.bp-9, 0x1.1p-10, -0x1.cp-10, 0x1.54p-9, -0x1p-12, -0x1.9cp-9,
+ 0x1.2p-10, -0x1.fp-10, 0x1.3p-9, -0x1.9p-11, 0x1.cp-9, 0x1p-12,
+ -0x1.88p-9, 0x1.28p-10, -0x1.2p-9, 0x1.fp-10, -0x1.9p-10, 0x1.4cp-9,
+ -0x1p-10, 0x1.9p-9, -0x1.2p-11, 0x1.c4p-9, -0x1p-12, 0x1.e8p-9,
+ -0x1p-14, 0x1.fcp-9, 0.0,
+};
+
+// Generated by Sollya with:
+// for i from 0 to 128 do {
+// r = 2^-8 * nearestint( 2^8 / (1 + i*2^-7) );
+// b = nearestint(log(r)*2^43) * 2^-43;
+// c = round(log(r) - b, D, RN);
+// print("{", -c, ",", -b, "},");
+// };
+// We replace LOG_R1_DD[128] with log(1.0) == 0.0
+constexpr fputil::DoubleDouble LOG_R1_DD[129] = {
+ {0.0, 0.0},
+ {-0x1.0c76b999d2be8p-46, 0x1.010157589p-7},
+ {-0x1.3dc5b06e2f7d2p-45, 0x1.0205658938p-6},
+ {-0x1.aa0ba325a0c34p-45, 0x1.8492528c9p-6},
+ {0x1.111c05cf1d753p-47, 0x1.0415d89e74p-5},
+ {-0x1.c167375bdfd28p-45, 0x1.466aed42ep-5},
+ {-0x1.29efbec19afa2p-47, 0x1.67c94f2d4cp-5},
+ {0x1.0fc1a353bb42ep-45, 0x1.aaef2d0fbp-5},
+ {-0x1.e113e4fc93b7bp-47, 0x1.eea31c006cp-5},
+ {-0x1.5325d560d9e9bp-45, 0x1.1973bd1466p-4},
+ {0x1.cc85ea5db4ed7p-45, 0x1.3bdf5a7d1ep-4},
+ {-0x1.53a2582f4e1efp-48, 0x1.4d3115d208p-4},
+ {0x1.c1e8da99ded32p-49, 0x1.700d30aeacp-4},
+ {0x1.3115c3abd47dap-45, 0x1.9335e5d594p-4},
+ {-0x1.e42b6b94407c8p-47, 0x1.a4e7640b1cp-4},
+ {0x1.646d1c65aacd3p-45, 0x1.c885801bc4p-4},
+ {0x1.a89401fa71733p-46, 0x1.da72763844p-4},
+ {-0x1.534d64fa10afdp-45, 0x1.fe89139dbep-4},
+ {0x1.1ef78ce2d07f2p-45, 0x1.1178e8227ep-3},
+ {0x1.ca78e44389934p-45, 0x1.1aa2b7e23fp-3},
+ {0x1.39d6ccb81b4a1p-47, 0x1.2d1610c868p-3},
+ {0x1.62fa8234b7289p-51, 0x1.365fcb0159p-3},
+ {0x1.5837954fdb678p-45, 0x1.4913d8333bp-3},
+ {0x1.633e8e5697dc7p-45, 0x1.527e5e4a1bp-3},
+ {-0x1.27023eb68981cp-46, 0x1.5bf406b544p-3},
+ {-0x1.5118de59c21e1p-45, 0x1.6f0128b757p-3},
+ {-0x1.c661070914305p-46, 0x1.7898d85445p-3},
+ {-0x1.73d54aae92cd1p-47, 0x1.8beafeb39p-3},
+ {0x1.7f22858a0ff6fp-47, 0x1.95a5adcf7p-3},
+ {0x1.9904d6865817ap-45, 0x1.9f6c407089p-3},
+ {-0x1.c358d4eace1aap-47, 0x1.b31d8575bdp-3},
+ {-0x1.d4bc4595412b6p-45, 0x1.bd087383bep-3},
+ {-0x1.1ec72c5962bd2p-48, 0x1.c6ffbc6f01p-3},
+ {-0x1.84a7e75b6f6e4p-47, 0x1.d1037f2656p-3},
+ {0x1.212276041f43p-51, 0x1.e530effe71p-3},
+ {-0x1.a211565bb8e11p-51, 0x1.ef5ade4ddp-3},
+ {0x1.bcbecca0cdf3p-46, 0x1.f991c6cb3bp-3},
+ {-0x1.6f08c1485e94ap-46, 0x1.01eae5626c8p-2},
+ {0x1.7188b163ceae9p-45, 0x1.0c42d67616p-2},
+ {-0x1.c210e63a5f01cp-45, 0x1.1178e8227e8p-2},
+ {0x1.b9acdf7a51681p-45, 0x1.16b5ccbacf8p-2},
+ {0x1.ca6ed5147bdb7p-45, 0x1.1bf99635a68p-2},
+ {0x1.a87deba46baeap-47, 0x1.214456d0eb8p-2},
+ {0x1.c93c1df5bb3b6p-45, 0x1.269621134d8p-2},
+ {0x1.a9cfa4a5004f4p-45, 0x1.2bef07cdc9p-2},
+ {0x1.16ecdb0f177c8p-46, 0x1.36b6776be1p-2},
+ {0x1.83b54b606bd5cp-46, 0x1.3c25277333p-2},
+ {0x1.8e436ec90e09dp-47, 0x1.419b423d5e8p-2},
+ {-0x1.f27ce0967d675p-45, 0x1.4718dc271c8p-2},
+ {-0x1.e20891b0ad8a4p-45, 0x1.4c9e09e173p-2},
+ {0x1.ebe708164c759p-45, 0x1.522ae0738ap-2},
+ {0x1.fadedee5d40efp-46, 0x1.57bf753c8dp-2},
+ {-0x1.a0b2a08a465dcp-47, 0x1.5d5bddf596p-2},
+ {-0x1.db623e731aep-45, 0x1.630030b3abp-2},
+ {0x1.0a0d32756ebap-45, 0x1.68ac83e9c68p-2},
+ {0x1.721657c222d87p-46, 0x1.6e60ee6af18p-2},
+ {0x1.d8b0949dc60b3p-45, 0x1.741d876c678p-2},
+ {0x1.9ec7d2efd1778p-45, 0x1.79e26687cf8p-2},
+ {-0x1.72090c812566ap-45, 0x1.7fafa3bd818p-2},
+ {0x1.fd56f3333778ap-45, 0x1.85855776dc8p-2},
+ {-0x1.05ae1e5e7047p-45, 0x1.8b639a88b3p-2},
+ {-0x1.766b52ee6307dp-46, 0x1.914a8635bf8p-2},
+ {-0x1.52313a502d9fp-46, 0x1.973a3431358p-2},
+ {-0x1.52313a502d9fp-46, 0x1.973a3431358p-2},
+ {-0x1.6279e10d0c0bp-45, 0x1.9d32bea15fp-2},
+ {0x1.3c6457f9d79f5p-45, 0x1.a33440224f8p-2},
+ {0x1.e36f2bea77a5dp-46, 0x1.a93ed3c8ad8p-2},
+ {-0x1.17cc552774458p-45, 0x1.af5295248dp-2},
+ {0x1.095252d841995p-46, 0x1.b56fa044628p-2},
+ {0x1.7d85bf40a666dp-45, 0x1.bb9611b80ep-2},
+ {0x1.cec807fe8e18p-45, 0x1.c1c60693fap-2},
+ {0x1.cec807fe8e18p-45, 0x1.c1c60693fap-2},
+ {-0x1.9b6ddc15249aep-45, 0x1.c7ff9c74558p-2},
+ {-0x1.797c33ec7a6bp-47, 0x1.ce42f180648p-2},
+ {0x1.35bafe9a767a8p-45, 0x1.d490246def8p-2},
+ {-0x1.ea42d60dc616ap-46, 0x1.dae75484c98p-2},
+ {-0x1.ea42d60dc616ap-46, 0x1.dae75484c98p-2},
+ {-0x1.326b207322938p-46, 0x1.e148a1a2728p-2},
+ {-0x1.465505372bd08p-45, 0x1.e7b42c3ddbp-2},
+ {0x1.f27f45a470251p-45, 0x1.ee2a156b41p-2},
+ {0x1.f27f45a470251p-45, 0x1.ee2a156b41p-2},
+ {0x1.2cde56f014a8bp-46, 0x1.f4aa7ee0318p-2},
+ {0x1.085fa3c164935p-47, 0x1.fb358af7a48p-2},
+ {-0x1.53ba3b1727b1cp-47, 0x1.00e5ae5b208p-1},
+ {-0x1.53ba3b1727b1cp-47, 0x1.00e5ae5b208p-1},
+ {-0x1.4c45fe79539ep-47, 0x1.04360be7604p-1},
+ {0x1.6812241edf5fdp-45, 0x1.078bf0533c4p-1},
+ {0x1.f486b887e7e27p-46, 0x1.0ae76e2d054p-1},
+ {0x1.f486b887e7e27p-46, 0x1.0ae76e2d054p-1},
+ {0x1.c299807801742p-46, 0x1.0e4898611ccp-1},
+ {-0x1.58647bb9ddcb2p-45, 0x1.11af823c75cp-1},
+ {-0x1.58647bb9ddcb2p-45, 0x1.11af823c75cp-1},
+ {-0x1.edd97a293ae49p-45, 0x1.151c3f6f298p-1},
+ {0x1.4cc4ef8ab465p-46, 0x1.188ee40f23cp-1},
+ {0x1.4cc4ef8ab465p-46, 0x1.188ee40f23cp-1},
+ {0x1.cacdeed70e667p-51, 0x1.1c07849ae6p-1},
+ {-0x1.a7242c9fe81d3p-45, 0x1.1f8635fc618p-1},
+ {-0x1.a7242c9fe81d3p-45, 0x1.1f8635fc618p-1},
+ {0x1.2fc066e48667bp-46, 0x1.230b0d8bebcp-1},
+ {-0x1.b61f10522625p-47, 0x1.269621134dcp-1},
+ {-0x1.b61f10522625p-47, 0x1.269621134dcp-1},
+ {0x1.06d2be797882dp-45, 0x1.2a2786d0ecp-1},
+ {-0x1.7a6e507b9dc11p-46, 0x1.2dbf557b0ep-1},
+ {-0x1.7a6e507b9dc11p-46, 0x1.2dbf557b0ep-1},
+ {-0x1.74e93c5a0ed9cp-45, 0x1.315da443408p-1},
+ {-0x1.74e93c5a0ed9cp-45, 0x1.315da443408p-1},
+ {0x1.0b83f9527e6acp-46, 0x1.35028ad9d8cp-1},
+ {-0x1.18b7abb5569a4p-45, 0x1.38ae2171978p-1},
+ {-0x1.18b7abb5569a4p-45, 0x1.38ae2171978p-1},
+ {-0x1.2b7367cfe13c2p-47, 0x1.3c6080c36cp-1},
+ {-0x1.2b7367cfe13c2p-47, 0x1.3c6080c36cp-1},
+ {-0x1.6ce7930f0c74cp-45, 0x1.4019c2125ccp-1},
+ {-0x1.6ce7930f0c74cp-45, 0x1.4019c2125ccp-1},
+ {-0x1.d984f481051f7p-48, 0x1.43d9ff2f924p-1},
+ {-0x1.2cb6af94d60aap-45, 0x1.47a1527e8a4p-1},
+ {-0x1.2cb6af94d60aap-45, 0x1.47a1527e8a4p-1},
+ {0x1.f7115ed4c541cp-49, 0x1.4b6fd6f970cp-1},
+ {0x1.f7115ed4c541cp-49, 0x1.4b6fd6f970cp-1},
+ {-0x1.e6c516d93b8fbp-45, 0x1.4f45a835a5p-1},
+ {-0x1.e6c516d93b8fbp-45, 0x1.4f45a835a5p-1},
+ {0x1.5ccc45d257531p-47, 0x1.5322e268678p-1},
+ {0x1.5ccc45d257531p-47, 0x1.5322e268678p-1},
+ {0x1.9980bff3303ddp-47, 0x1.5707a26bb8cp-1},
+ {0x1.9980bff3303ddp-47, 0x1.5707a26bb8cp-1},
+ {0x1.dfa63ac10c9fbp-45, 0x1.5af405c3648p-1},
+ {0x1.dfa63ac10c9fbp-45, 0x1.5af405c3648p-1},
+ {0x1.202380cda46bep-45, 0x1.5ee82aa2418p-1},
+ {0x1.202380cda46bep-45, 0x1.5ee82aa2418p-1},
+ {0.0, 0.0},
+};
+
+// Degree-7 minimax polynomial log(1 + v) ~ v - v^2 / 2 + ...
+// generated by Sollya with:
+// > P = fpminimax(log(1 + x)/x, 6, [|1, 1, D...|],
+// [-0x1.69000000000edp-8, 0x1.7f00000000081p-8]);
+constexpr double P_COEFFS[6] = {-0x1p-1,
+ 0x1.5555555555166p-2,
+ -0x1.fffffffdb7746p-3,
+ 0x1.99999a8718a6p-3,
+ -0x1.555874ce8ce22p-3,
+ 0x1.24335555ddbe5p-3};
+
+// -log(r1) with 128-bit precision generated by SageMath with:
+//
+// for i in range(129):
+// r = 2^-8 * round( 2^8 / (1 + i*2^(-7)) );
+// s, m, e = RealField(128)(r).log().sign_mantissa_exponent();
+// print("{false,", e, ", MType({", hex(m % 2^64), ",", hex((m >> 64) %
+// 2^64),
+// "})},");
+const Float128 LOG_R1[129] = {
+ {false, 0, MType(0)},
+ {false, -134, MType({0x662d417ced007a46, 0x8080abac46f38946})},
+ {false, -133, MType({0x91d082dce3ddcd38, 0x8102b2c49ac23a4f})},
+ {false, -133, MType({0xda5f3cc0b3251dbd, 0xc24929464655f45c})},
+ {false, -132, MType({0xb9e3aea6c444ef07, 0x820aec4f3a222380})},
+ {false, -132, MType({0x521016bd904dc968, 0xa33576a16f1f4c64})},
+ {false, -132, MType({0x27cca0bcc06c2f92, 0xb3e4a796a5dac208})},
+ {false, -132, MType({0xa9dda17056e45ed5, 0xd5779687d887e0d1})},
+ {false, -132, MType({0x606d89093278a939, 0xf7518e0035c3dd83})},
+ {false, -131, MType({0xa7c9859530a45153, 0x8cb9de8a32ab368a})},
+ {false, -131, MType({0x976d3b5b45f6ca0b, 0x9defad3e8f73217a})},
+ {false, -131, MType({0x3e858f08597b3a69, 0xa6988ae903f562ed})},
+ {false, -131, MType({0x6a677b4c8bec22e1, 0xb8069857560707a3})},
+ {false, -131, MType({0xeaf51f66692844ba, 0xc99af2eaca4c4570})},
+ {false, -131, MType({0x46bbf837b4d320c6, 0xd273b2058de1bd49})},
+ {false, -131, MType({0x196ab34ce0bccd12, 0xe442c00de2591b47})},
+ {false, -131, MType({0x3f4e2e660317d55f, 0xed393b1c22351280})},
+ {false, -131, MType({0xc17bd40d8d9291ec, 0xff4489cedeab2ca6})},
+ {false, -130, MType({0x9c5a0fe396f40f1e, 0x88bc74113f23def1})},
+ {false, -130, MType({0x88713268840cbcc0, 0x8d515bf11fb94f1c})},
+ {false, -130, MType({0x65c0da506a088484, 0x968b08643409ceb6})},
+ {false, -130, MType({0x411a5b944aca8708, 0x9b2fe580ac80b17d})},
+ {false, -130, MType({0xa9fb6cf0ecb411b7, 0xa489ec199dab06f2})},
+ {false, -130, MType({0xcad2fb8d48054ae0, 0xa93f2f250dac67d1})},
+ {false, -130, MType({0x149767e410316d2c, 0xadfa035aa1ed8fdc})},
+ {false, -130, MType({0x34c7bc3d32750fde, 0xb780945bab55dce4})},
+ {false, -130, MType({0x8f6ebcfb2016a439, 0xbc4c6c2a226399ef})},
+ {false, -130, MType({0xaa8b6997a402bf30, 0xc5f57f59c7f46155})},
+ {false, -130, MType({0x2c507fb7a3d0bf6a, 0xcad2d6e7b80bf914})},
+ {false, -130, MType({0xd0cb02f33f79c16c, 0xcfb6203844b3209a})},
+ {false, -130, MType({0x58a98f2ad65bee9b, 0xd98ec2bade71e539})},
+ {false, -130, MType({0x4d57da945b5d0aaa, 0xde8439c1dec56877})},
+ {false, -130, MType({0x4e9a750b6b68781d, 0xe37fde37807b84e3})},
+ {false, -130, MType({0xc524848e3443e040, 0xe881bf932af3dac0})},
+ {false, -130, MType({0x3b020fa1820c9492, 0xf29877ff38809091})},
+ {false, -130, MType({0x54d2238f75f969b1, 0xf7ad6f26e7ff2ef7})},
+ {false, -130, MType({0xca0cdf301431b60f, 0xfcc8e3659d9bcbec})},
+ {false, -129, MType({0xf5bd0b5b3479d5f4, 0x80f572b1363487b9})},
+ {false, -129, MType({0x163ceae88f720f1e, 0x86216b3b0b17188b})},
+ {false, -129, MType({0x9c5a0fe396f40f1e, 0x88bc74113f23def1})},
+ {false, -129, MType({0xf7a5168126a58b9a, 0x8b5ae65d67db9acd})},
+ {false, -129, MType({0x5147bdb6ddcaf59c, 0x8dfccb1ad35ca6ed})},
+ {false, -129, MType({0xae91aeba609c8877, 0x90a22b6875c6a1f7})},
+ {false, -129, MType({0xdf5bb3b60554e152, 0x934b1089a6dc93c1})},
+ {false, -129, MType({0x4a5004f3ef063313, 0x95f783e6e49a9cfa})},
+ {false, -129, MType({0xd878bbe3d392be25, 0x9b5b3bb5f088b766})},
+ {false, -129, MType({0x5b035eae273a855f, 0x9e1293b9998c1daa})},
+ {false, -129, MType({0xbb2438273918db7e, 0xa0cda11eaf46390d})},
+ {false, -129, MType({0xf698298adddd7f32, 0xa38c6e138e20d831})},
+ {false, -129, MType({0xe4f5275c2d15c21f, 0xa64f04f0b961df76})},
+ {false, -129, MType({0x8164c759686a2209, 0xa9157039c51ebe70})},
+ {false, -129, MType({0xf72ea07749ce6bd3, 0xabdfba9e468fd6f6})},
+ {false, -129, MType({0x7dd6e688ebb13b03, 0xaeadeefacaf97d35})},
+ {false, -129, MType({0x18ce51fff99479cd, 0xb1801859d56249dc})},
+ {false, -129, MType({0x2756eba00bc33978, 0xb45641f4e350a0d3})},
+ {false, -129, MType({0xbe1116c3466beb6d, 0xb730773578cb90b2})},
+ {false, -129, MType({0x49dc60b2b059a60b, 0xba0ec3b633dd8b09})},
+ {false, -129, MType({0x2efd17781bb3afec, 0xbcf13343e7d9ec7d})},
+ {false, -129, MType({0x37eda996244bccb0, 0xbfd7d1dec0a8df6f})},
+ {false, -129, MType({0x33337789d592e296, 0xc2c2abbb6e5fd56f})},
+ {false, -129, MType({0x1a18fb8f9f9ef280, 0xc5b1cd44596fa51e})},
+ {false, -129, MType({0x688ce7c1a75e341a, 0xc8a5431adfb44ca5})},
+ {false, -129, MType({0x2d7e9307c70c0668, 0xcb9d1a189ab56e76})},
+ {false, -129, MType({0x2d7e9307c70c0668, 0xcb9d1a189ab56e76})},
+ {false, -129, MType({0xef2f3f4f861ad6a9, 0xce995f50af69d861})},
+ {false, -129, MType({0x7f9d79f51dcc7301, 0xd19a201127d3c645})},
+ {false, -129, MType({0x5f53bd2e406e66e7, 0xd49f69e456cf1b79})},
+ {false, -129, MType({0xad88bba7d0cee8e0, 0xd7a94a92466e833a})},
+ {false, -129, MType({0x96c20cca6efe2ac5, 0xdab7d02231484a92})},
+ {false, -129, MType({0xf40a666c87842843, 0xddcb08dc0717d85b})},
+ {false, -129, MType({0x7fe8e1802aba24d6, 0xe0e30349fd1cec80})},
+ {false, -129, MType({0x7fe8e1802aba24d6, 0xe0e30349fd1cec80})},
+ {false, -129, MType({0x3eadb651b49ac53a, 0xe3ffce3a2aa64922})},
+ {false, -129, MType({0x304e1653e71d9973, 0xe72178c0323a1a0f})},
+ {false, -129, MType({0xe9a767a80d6d97e8, 0xea481236f7d35baf})},
+ {false, -129, MType({0x4f91cf4b33e42998, 0xed73aa4264b0ade9})},
+ {false, -129, MType({0x4f91cf4b33e42998, 0xed73aa4264b0ade9})},
+ {false, -129, MType({0xfc66eb6408ff6433, 0xf0a450d139366ca6})},
+ {false, -129, MType({0xac8d42f78d3e65d3, 0xf3da161eed6b9aaf})},
+ {false, -129, MType({0x5a470250d40ebe90, 0xf7150ab5a09f27f4})},
+ {false, -129, MType({0x5a470250d40ebe90, 0xf7150ab5a09f27f4})},
+ {false, -129, MType({0xb780a545a1b54dcf, 0xfa553f7018c966f2})},
+ {false, -129, MType({0x8f05924d258c14c5, 0xfd9ac57bd244217e})},
+ {false, -128, MType({0x89d1b09c70c4010a, 0x8072d72d903d588b})},
+ {false, -128, MType({0x89d1b09c70c4010a, 0x8072d72d903d588b})},
+ {false, -128, MType({0x30d58c3f7e2ea1f, 0x821b05f3b01d6774})},
+ {false, -128, MType({0x20f6fafe8fbb68b9, 0x83c5f8299e2b4091})},
+ {false, -128, MType({0xe21f9f89c1ab80b2, 0x8573b71682a7d21a})},
+ {false, -128, MType({0xe21f9f89c1ab80b2, 0x8573b71682a7d21a})},
+ {false, -128, MType({0x1e005d06dbfa8f8, 0x87244c308e670a66})},
+ {false, -128, MType({0x223111a707b6de2c, 0x88d7c11e3ad53cdc})},
+ {false, -128, MType({0x223111a707b6de2c, 0x88d7c11e3ad53cdc})},
+ {false, -128, MType({0x2eb628dba173c82d, 0x8a8e1fb794b09134})},
+ {false, -128, MType({0xbe2ad19415fe25a5, 0x8c47720791e53313})},
+ {false, -128, MType({0xbe2ad19415fe25a5, 0x8c47720791e53313})},
+ {false, -128, MType({0xbddae1ccce247838, 0x8e03c24d73003959})},
+ {false, -128, MType({0x9b00bf167e95da67, 0x8fc31afe30b2c6de})},
+ {false, -128, MType({0x9b00bf167e95da67, 0x8fc31afe30b2c6de})},
+ {false, -128, MType({0x9b92199ed1a4bab1, 0x918586c5f5e4bf01})},
+ {false, -128, MType({0xdf5bb3b60554e152, 0x934b1089a6dc93c1})},
+ {false, -128, MType({0xdf5bb3b60554e152, 0x934b1089a6dc93c1})},
+ {false, -128, MType({0xf3cbc416a2418012, 0x9513c36876083695})},
+ {false, -128, MType({0xbe1188fbc94e2f15, 0x96dfaabd86fa1646})},
+ {false, -128, MType({0xbe1188fbc94e2f15, 0x96dfaabd86fa1646})},
+ {false, -128, MType({0x1d2f89321647b358, 0x98aed221a03458b6})},
+ {false, -128, MType({0x1d2f89321647b358, 0x98aed221a03458b6})},
+ {false, -128, MType({0xe549f9aaea3cb5e1, 0x9a81456cec642e0f})},
+ {false, -128, MType({0xa2554b2dd4619e63, 0x9c5710b8cbb73a42})},
+ {false, -128, MType({0xa2554b2dd4619e63, 0x9c5710b8cbb73a42})},
+ {false, -128, MType({0x30603d87b6df81ad, 0x9e304061b5fda919})},
+ {false, -128, MType({0x30603d87b6df81ad, 0x9e304061b5fda919})},
+ {false, -128, MType({0x67879c5a30cd1242, 0xa00ce1092e5498c3})},
+ {false, -128, MType({0x67879c5a30cd1242, 0xa00ce1092e5498c3})},
+ {false, -128, MType({0xb7efae08e597e16, 0xa1ecff97c91e267b})},
+ {false, -128, MType({0x83594fab088c0d65, 0xa3d0a93f45169a4a})},
+ {false, -128, MType({0x83594fab088c0d65, 0xa3d0a93f45169a4a})},
+ {false, -128, MType({0xaf6a62a0dec6e073, 0xa5b7eb7cb860fb88})},
+ {false, -128, MType({0xaf6a62a0dec6e073, 0xa5b7eb7cb860fb88})},
+ {false, -128, MType({0x49362382a768847a, 0xa7a2d41ad270c9d7})},
+ {false, -128, MType({0x49362382a768847a, 0xa7a2d41ad270c9d7})},
+ {false, -128, MType({0x8ba4aea614d05701, 0xa991713433c2b998})},
+ {false, -128, MType({0x8ba4aea614d05701, 0xa991713433c2b998})},
+ {false, -128, MType({0x7fe6607ba902ef3c, 0xab83d135dc633301})},
+ {false, -128, MType({0x7fe6607ba902ef3c, 0xab83d135dc633301})},
+ {false, -128, MType({0xd60864fd949b4bd3, 0xad7a02e1b24efd31})},
+ {false, -128, MType({0xd60864fd949b4bd3, 0xad7a02e1b24efd31})},
+ {false, -128, MType({0x66d235ee63073dd, 0xaf74155120c9011c})},
+ {false, -128, MType({0x66d235ee63073dd, 0xaf74155120c9011c})},
+ {false, 0, MType(0)},
+};
+
+// Logarithm range reduction - Step 2:
+// s(k) = 2^-18 round( 2^18 / (1 + k*2^-14) ) - 1 for k = -91 .. 96
+// Output range:
+// [-0x1.1037c00000040271p-15 , 0x1.108480000008096cp-15]
+constexpr double S2[198] = {
+ 0x1.6ep-8, 0x1.6ap-8, 0x1.66p-8, 0x1.62p-8, 0x1.5dcp-8,
+ 0x1.59cp-8, 0x1.55cp-8, 0x1.51cp-8, 0x1.4dcp-8, 0x1.49cp-8,
+ 0x1.458p-8, 0x1.418p-8, 0x1.3d8p-8, 0x1.398p-8, 0x1.358p-8,
+ 0x1.318p-8, 0x1.2d8p-8, 0x1.294p-8, 0x1.254p-8, 0x1.214p-8,
+ 0x1.1d4p-8, 0x1.194p-8, 0x1.154p-8, 0x1.114p-8, 0x1.0dp-8,
+ 0x1.09p-8, 0x1.05p-8, 0x1.01p-8, 0x1.fap-9, 0x1.f2p-9,
+ 0x1.eap-9, 0x1.e2p-9, 0x1.d98p-9, 0x1.d18p-9, 0x1.c98p-9,
+ 0x1.c18p-9, 0x1.b98p-9, 0x1.b18p-9, 0x1.a98p-9, 0x1.a18p-9,
+ 0x1.998p-9, 0x1.91p-9, 0x1.89p-9, 0x1.81p-9, 0x1.79p-9,
+ 0x1.71p-9, 0x1.69p-9, 0x1.61p-9, 0x1.59p-9, 0x1.51p-9,
+ 0x1.49p-9, 0x1.41p-9, 0x1.388p-9, 0x1.308p-9, 0x1.288p-9,
+ 0x1.208p-9, 0x1.188p-9, 0x1.108p-9, 0x1.088p-9, 0x1.008p-9,
+ 0x1.f1p-10, 0x1.e1p-10, 0x1.d1p-10, 0x1.c1p-10, 0x1.b1p-10,
+ 0x1.a1p-10, 0x1.91p-10, 0x1.81p-10, 0x1.71p-10, 0x1.6p-10,
+ 0x1.5p-10, 0x1.4p-10, 0x1.3p-10, 0x1.2p-10, 0x1.1p-10,
+ 0x1p-10, 0x1.ep-11, 0x1.cp-11, 0x1.ap-11, 0x1.8p-11,
+ 0x1.6p-11, 0x1.4p-11, 0x1.2p-11, 0x1p-11, 0x1.cp-12,
+ 0x1.8p-12, 0x1.4p-12, 0x1p-12, 0x1.8p-13, 0x1p-13,
+ 0x1p-14, 0.0, -0x1p-14, -0x1p-13, -0x1.8p-13,
+ -0x1p-12, -0x1.4p-12, -0x1.8p-12, -0x1.cp-12, -0x1p-11,
+ -0x1.2p-11, -0x1.4p-11, -0x1.6p-11, -0x1.8p-11, -0x1.ap-11,
+ -0x1.cp-11, -0x1.ep-11, -0x1p-10, -0x1.1p-10, -0x1.2p-10,
+ -0x1.3p-10, -0x1.4p-10, -0x1.5p-10, -0x1.6p-10, -0x1.6fp-10,
+ -0x1.7fp-10, -0x1.8fp-10, -0x1.9fp-10, -0x1.afp-10, -0x1.bfp-10,
+ -0x1.cfp-10, -0x1.dfp-10, -0x1.efp-10, -0x1.ffp-10, -0x1.078p-9,
+ -0x1.0f8p-9, -0x1.178p-9, -0x1.1f8p-9, -0x1.278p-9, -0x1.2f8p-9,
+ -0x1.378p-9, -0x1.3fp-9, -0x1.47p-9, -0x1.4fp-9, -0x1.57p-9,
+ -0x1.5fp-9, -0x1.67p-9, -0x1.6fp-9, -0x1.77p-9, -0x1.7fp-9,
+ -0x1.87p-9, -0x1.8fp-9, -0x1.968p-9, -0x1.9e8p-9, -0x1.a68p-9,
+ -0x1.ae8p-9, -0x1.b68p-9, -0x1.be8p-9, -0x1.c68p-9, -0x1.ce8p-9,
+ -0x1.d68p-9, -0x1.dep-9, -0x1.e6p-9, -0x1.eep-9, -0x1.f6p-9,
+ -0x1.fep-9, -0x1.03p-8, -0x1.07p-8, -0x1.0bp-8, -0x1.0fp-8,
+ -0x1.12cp-8, -0x1.16cp-8, -0x1.1acp-8, -0x1.1ecp-8, -0x1.22cp-8,
+ -0x1.26cp-8, -0x1.2acp-8, -0x1.2e8p-8, -0x1.328p-8, -0x1.368p-8,
+ -0x1.3a8p-8, -0x1.3e8p-8, -0x1.428p-8, -0x1.464p-8, -0x1.4a4p-8,
+ -0x1.4e4p-8, -0x1.524p-8, -0x1.564p-8, -0x1.5a4p-8, -0x1.5ep-8,
+ -0x1.62p-8, -0x1.66p-8, -0x1.6ap-8, -0x1.6ep-8, -0x1.72p-8,
+ -0x1.75cp-8, -0x1.79cp-8, -0x1.7dcp-8,
+};
+
+// -log(r) for the second step, generated by SageMath with:
+//
+// for i in range(-91, 97):
+// r = 2^-18 * round( 2^18 / (1 + i*2^(-14)) );
+// s, m, e = RealField(128)(r).log().sign_mantissa_exponent();
+// print("{false," if (s == -1) else "{true,", e, ",
+// MType({", hex(m % 2^64), ",", hex((m >> 64) % 2^64), "})},");
+const Float128 LOG_R2[198] = {
+ {true, -135, MType({0xa0e061c5f7431c5e, 0xb67dab2a1a5742a4})},
+ {true, -135, MType({0x5d5bfe7b969ed6ec, 0xb4807f24af682939})},
+ {true, -135, MType({0x4d08702ddfabc23f, 0xb2834b35b4d54d5f})},
+ {true, -135, MType({0xd4d366508b9953df, 0xb0860f5ceba9be95})},
+ {true, -135, MType({0xac18a289f8f214a9, 0xae68f71aa09e8847})},
+ {true, -135, MType({0xd5b42054abb88c45, 0xac6baaeed676e8f1})},
+ {true, -135, MType({0x9809d58ee484964, 0xaa6e56d87cd632d6})},
+ {true, -135, MType({0xb9e6fc7c72f06d73, 0xa870fad754bb8791})},
+ {true, -135, MType({0x6f78d6d0105c00e2, 0xa67396eb1f231892})},
+ {true, -135, MType({0x28f712629209148, 0xa4762b139d0626e7})},
+ {true, -135, MType({0xc98d898ef172df02, 0xa258dfd10aedaa67})},
+ {true, -135, MType({0xfcc37c3c3062bfa1, 0xa05b63a373e60a83})},
+ {true, -135, MType({0x3eb450db05763c36, 0x9e5ddf89cf42f501})},
+ {true, -135, MType({0x7146a86fd458b775, 0x9c605383ddf1b88c})},
+ {true, -135, MType({0xc20a0c9281474436, 0x9a62bf9160dcb286})},
+ {true, -135, MType({0xcdc57316ec4aebc3, 0x986523b218eb4ed6})},
+ {true, -135, MType({0xc060dad74cef4273, 0x96677fe5c70207b9})},
+ {true, -135, MType({0xed8def1a3e433499, 0x9449f92d2ff44633})},
+ {true, -135, MType({0x3ce7a1f85c27b4fc, 0x924c45073220b5e0})},
+ {true, -135, MType({0xf2ca893449f7f2cb, 0x904e88f368fea63f})},
+ {true, -135, MType({0x8d77d9fabd2853cf, 0x8e50c4f1956699ed})},
+ {true, -135, MType({0x93e828d75b58ded4, 0x8c52f901782e20ec})},
+ {true, -135, MType({0x9f9605b053c5acf0, 0x8a552522d227d87a})},
+ {true, -135, MType({0x62a149393bca7241, 0x8857495564236ae0})},
+ {true, -135, MType({0xaea6b56ce89203d4, 0x86398719b66bac7c})},
+ {true, -135, MType({0x242bd86d00609b2, 0x843b9aef044e4dcc})},
+ {true, -135, MType({0xdaabf92774bac84e, 0x823da6d4c89c6927})},
+ {true, -135, MType({0xa1c6f3fc242ef8d0, 0x803faacac419abf2})},
+ {true, -136, MType({0xa225ebc02e6d9dd4, 0xfc834da16f0d9f57})},
+ {true, -136, MType({0xc33f6ad340ae18a9, 0xf88735ccc7433381})},
+ {true, -136, MType({0x70b2a4d38a242244, 0xf48b0e171249b6bc})},
+ {true, -136, MType({0x1d54819048b811b0, 0xf08ed67fd190e280})},
+ {true, -136, MType({0x9c21b650afe9ede0, 0xec52ca07ed95f236})},
+ {true, -136, MType({0x935519c96d30e463, 0xe85671adecd28aac})},
+ {true, -136, MType({0xba88f6f2e2672cfe, 0xe45a0970dc912ca7})},
+ {true, -136, MType({0xb1a8b84657ae069, 0xe05d91503e298bc8})},
+ {true, -136, MType({0xea3bff8d197b20a1, 0xdc61094b92ed70ef})},
+ {true, -136, MType({0xcdbb931d6fecc249, 0xd86471625c28b9e5})},
+ {true, -136, MType({0xd971d560d5f00820, 0xd467c9941b2158f5})},
+ {true, -136, MType({0x75563561244c090b, 0xd06b11e051175493})},
+ {true, -136, MType({0xdc393c9a3f3b380f, 0xcc6e4a467f44c6fa})},
+ {true, -136, MType({0xe6abe6e9e4ee2096, 0xc831a4c6f6fa709d})},
+ {true, -136, MType({0x3ce3c8228583a66e, 0xc434bc6124a0f16e})},
+ {true, -136, MType({0xb96a79f5c5a4963a, 0xc037c413c61bfd93})},
+ {true, -136, MType({0xaaef27337008679f, 0xbc3abbde5c8d9bde})},
+ {true, -136, MType({0xa49a3fcaddc8bc5a, 0xb83da3c06911e509})},
+ {true, -136, MType({0xe0254feb785362fa, 0xb4407bb96cbf035a})},
+ {true, -136, MType({0x9893a4e25ab9dc95, 0xb04343c8e8a53245})},
+ {true, -136, MType({0x5d8b0f40a3708915, 0xac45fbee5dcebe0b})},
+ {true, -136, MType({0x5f4c11c2c7a58c69, 0xa848a4294d40035d})},
+ {true, -136, MType({0xb348cc5df706ffba, 0xa44b3c7937f76efd})},
+ {true, -136, MType({0x9159f2c55a18befd, 0xa04dc4dd9eed7d60})},
+ {true, -136, MType({0xbdfdee41fe6a5a02, 0x9c1064563058bef3})},
+ {true, -136, MType({0x4580ddf89853254d, 0x9812cbe346475a24})},
+ {true, -136, MType({0xac75e10d61fc3ee8, 0x9415238353489ffb})},
+ {true, -136, MType({0xcad9b30b29736155, 0x90176b35d83ce8e2})},
+ {true, -136, MType({0x6f881deb98fc45f3, 0x8c19a2fa55fe9b14})},
+ {true, -136, MType({0x70a04b63b7248c96, 0x881bcad04d622a3e})},
+ {true, -136, MType({0xb4823fb48035eddd, 0x841de2b73f361722})},
+ {true, -136, MType({0x3364ccb5b13cd47f, 0x801feaaeac42ef38})},
+ {true, -137, MType({0xe306977b049f0ad5, 0xf843c56c2a969897})},
+ {true, -137, MType({0xe3c4d9e9619bc045, 0xf0479599f617a843})},
+ {true, -137, MType({0x4356d525b5e6432d, 0xe84b45e5bc76702c})},
+ {true, -137, MType({0x7839dcd7989339ab, 0xe04ed64e7f14697a})},
+ {true, -137, MType({0x4e21f045ecb76f23, 0xd85246d33f47230b})},
+ {true, -137, MType({0x902e248dd4ba9b28, 0xd0559772fe5840b0})},
+ {true, -137, MType({0xa44449067ef92e01, 0xc858c82cbd857a72})},
+ {true, -137, MType({0x17926207cc22e4e6, 0xc05bd8ff7e009bd2})},
+ {true, -137, MType({0x1c349622f3fa5d82, 0xb85ec9ea40ef8309})},
+ {true, -137, MType({0x97fa2fd0c9dc723e, 0xafe1c6ece1a058dd})},
+ {true, -137, MType({0x983e80897cf1e60f, 0xa7e47606048b1a65})},
+ {true, -137, MType({0x7199cd06ae5d39b3, 0x9fe705341d236102})},
+ {true, -137, MType({0x43cd18a72a051a96, 0x97e974762c5e8f58})},
+ {true, -137, MType({0x7b6d1248c3e1fd40, 0x8febc3cb332616ff})},
+ {true, -137, MType({0xf5572a8814c703af, 0x87edf332325777c5})},
+ {true, -138, MType({0x26828c92649a3a39, 0xffe0055455887de0})},
+ {true, -138, MType({0x82c550bd1216d82a, 0xefe3e4643a640cf3})},
+ {true, -138, MType({0xda6959f7f0e01bf0, 0xdfe7839214b4e8ae})},
+ {true, -138, MType({0xda93e2fa85a8f214, 0xcfeae2dbe5d6736d})},
+ {true, -138, MType({0xb47505bfa5a03b06, 0xbfee023faf0c2480})},
+ {true, -138, MType({0xb1475a5180a43520, 0xaff0e1bb718186ad})},
+ {true, -138, MType({0xa8740b91c95df537, 0x9ff3814d2e4a36b2})},
+ {true, -138, MType({0x57d895d35921b59c, 0x8ff5e0f2e661e1c6})},
+ {true, -139, MType({0x3c56c598c659c2a3, 0xfff0015535588833})},
+ {true, -139, MType({0x2ef8ec33ed9d782a, 0xdff3c0e497ea4eb1})},
+ {true, -139, MType({0x379eba7e6465ff63, 0xbff7008ff5e0c257})},
+ {true, -139, MType({0x3f972b783fcab757, 0x9ff9c0535073a370})},
+ {true, -140, MType({0xde026e271ee0549d, 0xfff8005551558885})},
+ {true, -140, MType({0xeceb47ea01f6c632, 0xbffb8023febc0c25})},
+ {true, -141, MType({0x7333c57857e1ed52, 0xfffc001554d55888})},
+ {true, -142, MType({0x87dde026fa704374, 0xfffe000555455588})},
+ {false, 0, MType({0x0, 0x0})},
+ {false, -141, MType({0x44999abe2fe2cc65, 0x80010002aab2aac4})},
+ {false, -140, MType({0x4eef381581464ccb, 0x8002000aaaeaac44})},
+ {false, -140, MType({0xdfeb485085f6f454, 0xc004802401440c26})},
+ {false, -139, MType({0x99abe3be3a1c6e93, 0x8004002aacaac445})},
+ {false, -139, MType({0x6bc1e20eac8448b4, 0xa00640535a37a37a})},
+ {false, -139, MType({0x979eedc064c242fd, 0xc00900900a20c275})},
+ {false, -139, MType({0xc72446cc1bf728bd, 0xe00c40e4bd6e4efd})},
+ {false, -138, MType({0xf381b821bbb569e5, 0x800800aabaac446e})},
+ {false, -138, MType({0x569b26aaa485ea5c, 0x900a20f319a3e273})},
+ {false, -138, MType({0x2dcf56c83c80b028, 0xa00c814d7c6a37f8})},
+ {false, -138, MType({0x5f69768284463b9b, 0xb00f21bbe3e388ee})},
+ {false, -138, MType({0xb48ea6c05e2773a1, 0xc0120240510c284c})},
+ {false, -138, MType({0x14d9d76196d8043a, 0xd01522dcc4f87991})},
+ {false, -138, MType({0xe016a611a4415d72, 0xe018839340d4f241})},
+ {false, -138, MType({0x661e135f49a47c40, 0xf01c2465c5e61b6f})},
+ {false, -137, MType({0xbe6bf0fa435e8383, 0x801002ab2ac4499a})},
+ {false, -137, MType({0x9a31ba0cbc030353, 0x881213337898871e})},
+ {false, -137, MType({0x54b57dfe0c4c840f, 0x901443cccd362c9f})},
+ {false, -137, MType({0x7ad1e9c315328f7e, 0x98169478296fad41})},
+ {false, -137, MType({0x1f3f686cf3d6be22, 0xa01905368e2389b3})},
+ {false, -137, MType({0xf105b66ec4703ede, 0xa81b9608fc3c50ec})},
+ {false, -137, MType({0x610848c68df4d233, 0xb01e46f074b0a0f3})},
+ {false, -137, MType({0x2e0efddf33a20464, 0xb7a0e9ed7613acb0})},
+ {false, -137, MType({0xc2cdb3c750f127b4, 0xbfa3d9008e042ffb})},
+ {false, -137, MType({0xbd9533786d3f4c49, 0xc7a6e82ba36a7073})},
+ {false, -137, MType({0x82e237c9a4d450e3, 0xcfaa176fb76c8eb1})},
+ {false, -137, MType({0xc00b46a4d0e3dfd0, 0xd7ad66cdcb3cbe14})},
+ {false, -137, MType({0xea999c0df8546710, 0xdfb0d646e0194584})},
+ {false, -137, MType({0xcec6c2a9ad974f4f, 0xe7b465dbf74c8032})},
+ {false, -137, MType({0x2d2045da1570a07c, 0xefb8158e122cde5a})},
+ {false, -137, MType({0x6752e9b2381e3edc, 0xf7bbe55e321ce603})},
+ {false, -137, MType({0x3c1ed52728e00e40, 0xffbfd54d588b33c5})},
+ {false, -136, MType({0x493b0d873fb9a340, 0x83e1f2ae43793dc3})},
+ {false, -136, MType({0x29e38750c9d26893, 0x87e40ac65f6cc4a0})},
+ {false, -136, MType({0xaab9e8327258ac3f, 0x8be632ef80e9a0df})},
+ {false, -136, MType({0x28bc403d8a5f3c63, 0x8fe86b2a28bf51b3})},
+ {false, -136, MType({0xf720c1c97227fcdc, 0x93eab376d7c36377})},
+ {false, -136, MType({0x6ad9a3e3d11b66c1, 0x97ed0bd60ed17018})},
+ {false, -136, MType({0xedb27b79c90b4019, 0x9bef74484ecb1f6c})},
+ {false, -136, MType({0xa092a0d7ab21722a, 0x9fb1c4cd27012e19})},
+ {false, -136, MType({0x535d52f0939a4d02, 0xa3b44c65b71c2d85})},
+ {false, -136, MType({0x90a57e11edc1864e, 0xa7b6e412cadcb3dc})},
+ {false, -136, MType({0x68e9c90160031159, 0xabb98bd4e33c4381})},
+ {false, -136, MType({0xbf60594f929adeb8, 0xafbc43ac813a6ea3})},
+ {false, -136, MType({0x8a42158886775205, 0xb3bf0b9a25dcd7a2})},
+ {false, -136, MType({0x1ab45417663dee9e, 0xb7c1e39e522f316d})},
+ {false, -136, MType({0x6c51ae3ce1aea68a, 0xbbc4cbb987433fe4})},
+ {false, -136, MType({0x7c52ae8b40ebabb7, 0xbfc7c3ec4630d83c})},
+ {false, -136, MType({0xa857126f7cfaaa67, 0xc3cacc371015e15d})},
+ {false, -136, MType({0x14d05662cd29464a, 0xc7cde49a66165446})},
+ {false, -136, MType({0x8379db06ef3cd6bb, 0xcb90da1644d29bb7})},
+ {false, -136, MType({0x9025f4c67dd38bb6, 0xcf9411aa99ddb7de})},
+ {false, -136, MType({0xd6f8a61c892032ee, 0xd3975958f681086d})},
+ {false, -136, MType({0x9a2f20b4e2332d47, 0xd79ab121dbf8714c})},
+ {false, -136, MType({0x3c767d61f51d375b, 0xdb9e1905cb85ea59})},
+ {false, -136, MType({0xd4b2bd65bb25493c, 0xdfa1910546717fca})},
+ {false, -136, MType({0xc96c1254a30ef91f, 0xe3a51920ce095292})},
+ {false, -136, MType({0x73e324ce0946b214, 0xe7a8b158e3a198be})},
+ {false, -136, MType({0xcacd125a12bac62c, 0xebac59ae08949dd8})},
+ {false, -136, MType({0xcafdc27227b71eaa, 0xef6fd620b2b7a503})},
+ {false, -136, MType({0x688d4282f6026aa3, 0xf3739daf959aaafc})},
+ {false, -136, MType({0xe54e9e3804464cdd, 0xf777755d03f4e0b6})},
+ {false, -136, MType({0xcb78b383f4b59dce, 0xfb7b5d297f388a12})},
+ {false, -136, MType({0xee055fc515062c04, 0xff7f551588de024f})},
+ {false, -135, MType({0x207812b43382acdd, 0x81c1ae90d131de38})},
+ {false, -135, MType({0xdc90c4c4b61f3a87, 0x83c3baa726a721cc})},
+ {false, -135, MType({0x1a03f13fb2c978b1, 0x85c5cece05941dbc})},
+ {false, -135, MType({0xb36f282e83a7dc36, 0x87c7eb05aec1304f})},
+ {false, -135, MType({0xd82a46616d4c393f, 0x89a9eccd56a980c0})},
+ {false, -135, MType({0xbc6ff84713c9babd, 0x8bac18a640185360})},
+ {false, -135, MType({0x9f7942a516fc2d8a, 0x8dae4c90b22574f4})},
+ {false, -135, MType({0x15e50cfd9b29b427, 0x8fb0888ceda546ab})},
+ {false, -135, MType({0x9f465296ae7dd49a, 0x91b2cc9b336f3718})},
+ {false, -135, MType({0xb49c1eb9b348e6e4, 0x93b518bbc45dc268})},
+ {false, -135, MType({0xdaa320cd64c9d9c7, 0x95b76ceee14e728e})},
+ {false, -135, MType({0x75a91950ffe1e3b5, 0x9799a333de49b963})},
+ {false, -135, MType({0x5c6abcbf43f03f14, 0x999c070ba32068cd})},
+ {false, -135, MType({0x5a9e7f265d1ed157, 0x9b9e72f6b295ad4f})},
+ {false, -135, MType({0xefeb98d02a195c17, 0x9da0e6f54d9318fd})},
+ {false, -135, MType({0x2aa503a3110ab5a7, 0x9fa36307b5054ca8})},
+ {false, -135, MType({0xd0fe7e05869eb825, 0xa1a5e72e29dbf808})},
+ {false, -135, MType({0xe80a28f4e1e500d2, 0xa3884a68a750cb10})},
+ {false, -135, MType({0x531064151ca6e30b, 0xa58ade36aeef9f0b})},
+ {false, -135, MType({0x27c01ffa8e2e3c4b, 0xa78d7a1982c4b08f})},
+ {false, -135, MType({0x7ba9408dc857d568, 0xa9901e1163cbbbf5})},
+ {false, -135, MType({0x104d1e3331d3b4fa, 0xab92ca1e93038d76})},
+ {false, -135, MType({0x9343c846fcdf9137, 0xad957e41516e0158})},
+ {false, -135, MType({0x3977e89aec59bfa2, 0xaf780e79b2514889})},
+ {false, -135, MType({0x913d4e3dc55c3e6e, 0xb17ad246ef3713bc})},
+ {false, -135, MType({0x777b52a9e70d8bcc, 0xb37d9e2a7a56b09d})},
+ {false, -135, MType({0x55de916fd30591de, 0xb580722494be0c91})},
+ {false, -135, MType({0xe79cfb37be2861e4, 0xb7834e357f7e2600})},
+ {false, -135, MType({0x90983104d3805389, 0xb986325d7bab0c89})},
+ {false, -135, MType({0x59e3b2ec71ce64f4, 0xbb68ef9c254aa378})},
+ {false, -135, MType({0xe83183bf3dd612ef, 0xbd6be3718c77636f})},
+ {false, -135, MType({0xc4e3b0ac2fd52b7f, 0xbf6edf5ec44d9d35})},
+};
+
+// Logarithm range reduction - Step 3:
+// s(k) = 2^-21 round( 2^21 / (1 + k*2^-21) ) - 1 for k = -69 .. 69
+// Output range:
+// [-0x1.012bb800000800114p-22, 0x1p-22 ]
+constexpr double S3[139] = {
+ 0x1.14p-15, 0x1.1p-15, 0x1.0cp-15, 0x1.08p-15, 0x1.04p-15, 0x1p-15,
+ 0x1.f8p-16, 0x1.fp-16, 0x1.e8p-16, 0x1.ep-16, 0x1.d8p-16, 0x1.dp-16,
+ 0x1.c8p-16, 0x1.cp-16, 0x1.b8p-16, 0x1.bp-16, 0x1.a8p-16, 0x1.ap-16,
+ 0x1.98p-16, 0x1.9p-16, 0x1.88p-16, 0x1.8p-16, 0x1.78p-16, 0x1.7p-16,
+ 0x1.68p-16, 0x1.6p-16, 0x1.58p-16, 0x1.5p-16, 0x1.48p-16, 0x1.4p-16,
+ 0x1.38p-16, 0x1.3p-16, 0x1.28p-16, 0x1.2p-16, 0x1.18p-16, 0x1.1p-16,
+ 0x1.08p-16, 0x1p-16, 0x1.fp-17, 0x1.ep-17, 0x1.dp-17, 0x1.cp-17,
+ 0x1.bp-17, 0x1.ap-17, 0x1.9p-17, 0x1.8p-17, 0x1.7p-17, 0x1.6p-17,
+ 0x1.5p-17, 0x1.4p-17, 0x1.3p-17, 0x1.2p-17, 0x1.1p-17, 0x1p-17,
+ 0x1.ep-18, 0x1.cp-18, 0x1.ap-18, 0x1.8p-18, 0x1.6p-18, 0x1.4p-18,
+ 0x1.2p-18, 0x1p-18, 0x1.cp-19, 0x1.8p-19, 0x1.4p-19, 0x1p-19,
+ 0x1.8p-20, 0x1p-20, 0x1p-21, 0.0, -0x1p-21, -0x1p-20,
+ -0x1.8p-20, -0x1p-19, -0x1.4p-19, -0x1.8p-19, -0x1.cp-19, -0x1p-18,
+ -0x1.2p-18, -0x1.4p-18, -0x1.6p-18, -0x1.8p-18, -0x1.ap-18, -0x1.cp-18,
+ -0x1.ep-18, -0x1p-17, -0x1.1p-17, -0x1.2p-17, -0x1.3p-17, -0x1.4p-17,
+ -0x1.5p-17, -0x1.6p-17, -0x1.7p-17, -0x1.8p-17, -0x1.9p-17, -0x1.ap-17,
+ -0x1.bp-17, -0x1.cp-17, -0x1.dp-17, -0x1.ep-17, -0x1.fp-17, -0x1p-16,
+ -0x1.08p-16, -0x1.1p-16, -0x1.18p-16, -0x1.2p-16, -0x1.28p-16, -0x1.3p-16,
+ -0x1.38p-16, -0x1.4p-16, -0x1.48p-16, -0x1.5p-16, -0x1.58p-16, -0x1.6p-16,
+ -0x1.68p-16, -0x1.7p-16, -0x1.78p-16, -0x1.8p-16, -0x1.88p-16, -0x1.9p-16,
+ -0x1.98p-16, -0x1.ap-16, -0x1.a8p-16, -0x1.bp-16, -0x1.b8p-16, -0x1.cp-16,
+ -0x1.c8p-16, -0x1.dp-16, -0x1.d8p-16, -0x1.ep-16, -0x1.e8p-16, -0x1.fp-16,
+ -0x1.f8p-16, -0x1p-15, -0x1.04p-15, -0x1.08p-15, -0x1.0cp-15, -0x1.1p-15,
+ -0x1.14p-15,
+};
+
+// -log(r) for the third step, generated by SageMath with:
+//
+// for i in range(-69, 70):
+// r = 2^-21 * round( 2^21 / (1 + i*2^(-21)) );
+// s, m, e = RealField(128)(r).log().sign_mantissa_exponent();
+// print("{false," if (s == -1) else "{true,", e, ",
+// MType({", hex(m % 2^64), ",", hex((m >> 64) % 2^64), "})},");
+const Float128 LOG_R3[139] = {
+ {true, -142, MType({0xe39d3faf42340ed7, 0x89ff6b38d5de2622})},
+ {true, -142, MType({0x7ff3326682c02485, 0x87ff6f80ccb40f16})},
+ {true, -142, MType({0x5caf4fbe343cf928, 0x85ff73b8c3cdf731})},
+ {true, -142, MType({0xcdb6e554348f7fe8, 0x83ff77e0bb2ade79})},
+ {true, -142, MType({0xef009c2457de25d, 0x81ff7bf8b2c9c4f6})},
+ {true, -143, MType({0x8883333c57b57c74, 0xffff000155535558})},
+ {true, -143, MType({0xf32668f39c70d183, 0xfbff07f145931f44})},
+ {true, -143, MType({0x459a73c6a6486fe3, 0xf7ff0fc13650e7bd})},
+ {true, -143, MType({0x37b18cca7dd3a29f, 0xf3ff1771278aaecd})},
+ {true, -143, MType({0x513f610d21bcfc78, 0xefff1f01193e7480})},
+ {true, -143, MType({0xea190b95c0690b7b, 0xebff26710b6a38e1})},
+ {true, -143, MType({0x2a150f64f0ad1743, 0xe7ff2dc0fe0bfbfd})},
+ {true, -143, MType({0x90b5174e995e9d1, 0xe3ff34f0f121bddd})},
+ {true, -143, MType({0x4ed512b9b93ea2bf, 0xdfff3c00e4a97e8c})},
+ {true, -143, MType({0x934cea217ab794a2, 0xdbff42f0d8a13e15})},
+ {true, -143, MType({0x3e4ebe948afd2c76, 0xd7ff49c0cd06fc83})},
+ {true, -143, MType({0x87b7c0f5bcfee2e1, 0xd3ff5070c1d8b9df})},
+ {true, -143, MType({0x776666228cb6371b, 0xcfff5700b7147634})},
+ {true, -143, MType({0xe53a60f3514db358, 0xcbff5d70acb8318b})},
+ {true, -143, MType({0x79149c3b6e57fa86, 0xc7ff63c0a2c1ebef})},
+ {true, -143, MType({0xaad734c98416df2a, 0xc3ff69f0992fa568})},
+ {true, -143, MType({0xc26573679ed28334, 0xbfff70008fff5e00})},
+ {true, -143, MType({0xd7a3c6db6540809f, 0xbbff75f0872f15c0})},
+ {true, -143, MType({0xd277bde645fb1aad, 0xb7ff7bc07ebcccb1})},
+ {true, -143, MType({0x6ac80145a4087793, 0xb3ff817076a682dc})},
+ {true, -143, MType({0x287c4db30271e265, 0xafff87006eea3849})},
+ {true, -143, MType({0x637d6de42eeb151e, 0xabff8c706785ed00})},
+ {true, -143, MType({0x43b5348b6b898a8c, 0xa7ff91c06077a10a})},
+ {true, -143, MType({0xc10e7657978bd7f6, 0xa3ff96f059bd546e})},
+ {true, -143, MType({0xa37503f457310e59, 0x9fff9c0053550735})},
+ {true, -143, MType({0x82d5a40a3aa022ff, 0x9bffa0f04d3cb966})},
+ {true, -143, MType({0xc71e0d3ee3df5f4d, 0x97ffa5c047726b08})},
+ {true, -143, MType({0xa83ce0352bdbd79b, 0x93ffaa7041f41c23})},
+ {true, -143, MType({0x2e21a18d4680e8e4, 0x8fffaf003cbfccbe})},
+ {true, -143, MType({0x30bcb3e4e5dfbd28, 0x8bffb37037d37cdf})},
+ {true, -143, MType({0x57ff51d75c66d64a, 0x87ffb7c0332d2c8d})},
+ {true, -143, MType({0x1bdb87fdbe299f43, 0x83ffbbf02ecadbcf})},
+ {true, -144, MType({0x88885dde02700703, 0xffff800055551555})},
+ {true, -144, MType({0xd259ca803a0c1870, 0xf7ff87e04d94724c})},
+ {true, -144, MType({0xe514130851c7070a, 0xefff8f80464fce8f})},
+ {true, -144, MType({0x30a16898f3073a64, 0xe7ff96e03f832a2a})},
+ {true, -144, MType({0xc4ed64517b2949ce, 0xdfff9e00392a8526})},
+ {true, -144, MType({0x51e4fb4e32cf6350, 0xd7ffa4e03341df90})},
+ {true, -144, MType({0x277672a88350bcce, 0xcfffab802dc53971})},
+ {true, -144, MType({0x359153772a490f06, 0xc7ffb1e028b092d3})},
+ {true, -144, MType({0xc265ece6b481a0e, 0xbfffb80023ffebc0})},
+ {true, -144, MType({0xdb2781c03fa132f6, 0xb7ffbde01faf4440})},
+ {true, -144, MType({0x7287c95c845ada33, 0xafffc3801bba9c5e})},
+ {true, -144, MType({0x423b56b1263e5a77, 0xa7ffc8e0181df421})},
+ {true, -144, MType({0x5a3752ca4c076fa3, 0x9fffce0014d54b91})},
+ {true, -144, MType({0x6a71e2b27eb3f573, 0x97ffd2e011dca2b6})},
+ {true, -144, MType({0xc2e21b72cff39d8f, 0x8fffd7800f2ff997})},
+ {true, -144, MType({0x537ff612feb7ac9e, 0x87ffdbe00ccb503c})},
+ {true, -145, MType({0x5888873333c57c18, 0xffffc00015554d55})},
+ {true, -145, MType({0xfa51421842311c42, 0xefffc7c01193f9d1})},
+ {true, -145, MType({0x2c4ed6de475b942c, 0xdfffcf000e4aa5fa})},
+ {true, -145, MType({0xce77678cbb6fcb88, 0xcfffd5c00b7151d8})},
+ {true, -145, MType({0xc26629a679ed3b, 0xbfffdc0008fffd78})},
+ {true, -145, MType({0x23287cb9d3072728, 0xafffe1c006eea8e1})},
+ {true, -145, MType({0xd5a37540fd057315, 0x9fffe7000535541c})},
+ {true, -145, MType({0xf82e21c1fce36810, 0x8fffebc003cbff32})},
+ {true, -146, MType({0x5588887ddde02702, 0xffffe00005555455})},
+ {true, -146, MType({0x9ac4ed72adf5b295, 0xdfffe7800392aa14})},
+ {true, -146, MType({0xc26648066b482, 0xbfffee00023fffaf})},
+ {true, -146, MType({0x455a3754b292c077, 0x9ffff380014d552e})},
+ {true, -147, MType({0x5558888833333c58, 0xfffff00001555535})},
+ {true, -147, MType({0xe000c2665736679f, 0xbffff700008ffff5})},
+ {true, -148, MType({0x5555888885ddde02, 0xfffff80000555551})},
+ {true, -149, MType({0xd555588888733334, 0xfffffc0000155554})},
+ {false, 0, MType({0x0, 0x0})},
+ {false, -148, MType({0xeaaaac44444eeeef, 0x80000200000aaaaa})},
+ {false, -147, MType({0xaaaac444459999ac, 0x80000400002aaaac})},
+ {false, -147, MType({0x2000c2667596679f, 0xc00009000090000a})},
+ {false, -146, MType({0xaaac44446eeef381, 0x8000080000aaaaba})},
+ {false, -146, MType({0x655a3755f81815cc, 0xa0000c80014d557c})},
+ {false, -146, MType({0xc26684c66b482, 0xc000120002400051})},
+ {false, -146, MType({0xbac4ed7c40fb07eb, 0xe00018800392ab40})},
+ {false, -145, MType({0xaac44449999abe2c, 0x8000100002aaab2a})},
+ {false, -145, MType({0x82e21d79cbb6812, 0x9000144003cc00cd})},
+ {false, -145, MType({0xd5a37569adb01dc3, 0xa00019000535568d})},
+ {false, -145, MType({0x33287d01e8c9d1d9, 0xb0001e4006eeac74})},
+ {false, -145, MType({0xc266a32679ed48, 0xc000240009000288})},
+ {false, -145, MType({0xde77685122b2764b, 0xd0002a400b7158d1})},
+ {false, -145, MType({0x2c4ed810a8063f03, 0xe00031000e4aaf5b})},
+ {false, -145, MType({0xa5143e7be891c8f, 0xf00038401194062e})},
+ {false, -144, MType({0xac4444eeef3813a1, 0x800020000aaaaeaa})},
+ {false, -144, MType({0x5b7ff7fe1339025b, 0x880024200ccb5a6e})},
+ {false, -144, MType({0x42e21e26caf39e33, 0x900028800f300668})},
+ {false, -144, MType({0xf271e66fa5554bc6, 0x98002d2011dcb29e})},
+ {false, -144, MType({0x5a3757e0615cc676, 0xa000320014d55f19})},
+ {false, -144, MType({0xca3b5d8210ca5cab, 0xa8003720181e0bde})},
+ {false, -144, MType({0xf287d25f3cb032bb, 0xb0003c801bbab8f6})},
+ {false, -144, MType({0xe3278d840be28cdb, 0xb80042201faf6669})},
+ {false, -144, MType({0xc266dfe6b482076, 0xc000480024001440})},
+ {false, -144, MType({0x3d9166de380a6d3d, 0xc8004e2028b0c282})},
+ {false, -144, MType({0xa7768b356ba61e4b, 0xd00054802dc57139})},
+ {false, -144, MType({0xd9e51a1849db73c1, 0xd8005b203342206f})},
+ {false, -144, MType({0xc4ed8a9d907eb521, 0xe0006200392ad02e})},
+ {false, -144, MType({0xb8a197dea928acd7, 0xe80069203f838080})},
+ {false, -144, MType({0x65144cf7dcc72d3b, 0xf000708046503170})},
+ {false, -144, MType({0xda5a1108890d9f6a, 0xf80078204d94e308})},
+ {false, -143, MType({0xc4445999abe2ce2c, 0x800040002aaacaaa})},
+ {false, -143, MType({0x1fdbbb4f3bffc832, 0x840044102ecb2431})},
+ {false, -143, MType({0x97ff8f39ec91b4ee, 0x88004840332d7e1d})},
+ {false, -143, MType({0x74bcfcf0b3f0a95d, 0x8c004c9037d3d876})},
+ {false, -143, MType({0x2e21f80ca6813aff, 0x900051003cc03342})},
+ {false, -143, MType({0x6c3d4629170ce87f, 0x9400559041f48e87})},
+ {false, -143, MType({0x71e84e3b80a8881, 0x98005a404772ea4d})},
+ {false, -143, MType({0x6d62fdcbdd6bec3, 0x9c005f104d3d469a})},
+ {false, -143, MType({0xa375a6b701dc77c0, 0xa00064005355a375})},
+ {false, -143, MType({0x450f331826ad6b05, 0xa400691059be00e7})},
+ {false, -143, MType({0x83b60ea8bd0aa459, 0xa8006e4060785ef6})},
+ {false, -143, MType({0x277e691469dd13f5, 0xac0073906786bdab})},
+ {false, -143, MType({0x287d6e0a0d1e25eb, 0xb00079006eeb1d0d})},
+ {false, -143, MType({0xaec94b3be9b060f5, 0xb4007e9076a77d24})},
+ {false, -143, MType({0x1279365fce280cce, 0xb80084407ebdddfa})},
+ {false, -143, MType({0xdba5732f3e83e04a, 0xbc008a1087303f95})},
+ {false, -143, MType({0xc26759679ed5b754, 0xc00090009000a200})},
+ {false, -143, MType({0xaed95aca5edb5109, 0xc400961099310543})},
+ {false, -143, MType({0xb917091d2687160f, 0xc8009c40a2c36967})},
+ {false, -143, MType({0x293d1c2a0378e75d, 0xcc00a290acb9ce76})},
+ {false, -143, MType({0x776977bf9766f5a7, 0xd000a900b7163478})},
+ {false, -143, MType({0x4bbb31b14776a18b, 0xd400af90c1da9b78})},
+ {false, -143, MType({0x7e5297d76c8564ba, 0xd800b640cd09037f})},
+ {false, -143, MType({0x1751360f8461c447, 0xdc00bd10d8a36c98})},
+ {false, -143, MType({0x4ed9dc3c63f44c41, 0xe000c400e4abd6cc})},
+ {false, -143, MType({0x8d10a4466a5894d5, 0xe400cb10f1244226})},
+ {false, -143, MType({0x6a1af81bb4e6510e, 0xe800d240fe0eaeb1})},
+ {false, -143, MType({0xae1f97b0542a677a, 0xec00d9910b6d1c77})},
+ {false, -143, MType({0x51469efe81d014cc, 0xf000e10119418b84})},
+ {false, -143, MType({0x7bb98c06d77a18b4, 0xf400e891278dfbe2})},
+ {false, -143, MType({0x85a344d0868bed17, 0xf800f04136546d9d})},
+ {false, -143, MType({0xf7301d6990e307cc, 0xfc00f8114596e0c0})},
+ {false, -142, MType({0x4446eef38140138f, 0x80008000aaabaaac})},
+ {false, -142, MType({0x10f5e43296105497, 0x82008408b2cbe5b8})},
+ {false, -142, MType({0xedbd4f83ef63f730, 0x84008820bb2d2189})},
+ {false, -142, MType({0xfeb654fd541c638e, 0x86008c48c3d05e27})},
+ {false, -142, MType({0x7ffadeb8882f7674, 0x88009080ccb69b98})},
+ {false, -142, MType({0xc5a59fd36bd44397, 0x8a0094c8d5e0d9e1})},
+};
+
+// Minimax polynomial generated by Sollya with:
+// > P = fpminimax((log(1 + x) - x)/x^2, 3, [|1, 128...|],
+// [-0x1.01928p-22 , 0x1p-22]);
+// > P;
+// > dirtyinfnorm(log(1 + x)/x - 1 - x*P, [-0x1.01928p-22 , 0x1p-22]);
+// 0x1.ce1e...p-116
+const Float128 BIG_COEFFS[4]{
+ {false, -130, MType({0x7ed78465d460315b, 0xccccccd74818e397})},
+ {true, -129, MType({0xc6388a23871ce156, 0x80000000000478b0})},
+ {false, -129, MType({0xaa807bd867763262, 0xaaaaaaaaaaaaaaaa})},
+ {true, -128, MType({0x0, 0x8000000000000000})},
+};
+
+LIBC_INLINE double log1p_accurate(int e_x, int index,
+ fputil::DoubleDouble m_x) {
+ Float128 e_x_f128(static_cast<float>(e_x));
+ Float128 sum = fputil::quick_mul(LOG_2, e_x_f128);
+ sum = fputil::quick_add(sum, LOG_R1[index]);
+
+ // fputil::DoubleDouble v4;
+ Float128 v = fputil::quick_add(Float128(m_x.hi), Float128(m_x.lo));
+
+ // Skip 2nd range reduction step if |m_x| <= 2^-15.
+ if (m_x.hi > 0x1p-15 || m_x.hi < -0x1p-15) {
+ // Range reduction - Step 2.
+ // For k such that: k * 2^-14 - 2^-15 <= m_x.hi < k * 2^-14 + 2^-15,
+ // Let s_k = 2^-18 * round( 2^18 / (1 + k*2^-14) ) - 1
+ // Then the 2nd reduced argument is:
+ // (1 + s_k) * (1 + m_x) - 1 =
+ // = s_k + m_x + s_k * m_x
+ // Output range:
+ // -0x1.1037c00000040271p-15 <= v2.hi + v2.lo <= 0x1.108480000008096cp-15
+ int idx2 = static_cast<int>(0x1p14 * (m_x.hi + (91 * 0x1p-14 + 0x1p-15)));
+ sum = fputil::quick_add(sum, LOG_R2[idx2]);
+ Float128 s2 = Float128(S2[idx2]);
+ v = fputil::quick_add(fputil::quick_add(v, s2), fputil::quick_mul(v, s2));
+ }
+
+ // Skip 3rd range reduction step if |v| <= 2^-22.
+ if (v.exponent > -150) {
+ // Range reduction - Step 3.
+ // For k such that: k * 2^-21 - 2^-22 <= v2.hi < k * 2^-21 + 2^-22,
+ // Let s_k = 2^-21 * round( 2^21 / (1 + k*2^-21) ) - 1
+ // Then the 3rd reduced argument is:
+ // v3.hi + v3.lo ~ (1 + s_k) * (1 + v2.hi + v2.lo) - 1
+ // Output range:
+ // -0x1.012bb800000800114p-22 <= v3.hi + v3.lo <= 0x1p-22
+ int idx3 =
+ static_cast<int>(0x1p21 * (double(v) + (69 * 0x1p-21 + 0x1p-22)));
+ sum = fputil::quick_add(sum, LOG_R3[idx3]);
+ Float128 s3 = Float128(S3[idx3]);
+ v = fputil::quick_add(fputil::quick_add(v, s3), fputil::quick_mul(v, s3));
+ }
+
+ // Polynomial approximation
+ Float128 p = fputil::quick_mul(v, BIG_COEFFS[0]);
+ p = fputil::quick_mul(v, fputil::quick_add(p, BIG_COEFFS[1]));
+ p = fputil::quick_mul(v, fputil::quick_add(p, BIG_COEFFS[2]));
+ p = fputil::quick_mul(v, fputil::quick_add(p, BIG_COEFFS[3]));
+ p = fputil::quick_add(v, fputil::quick_mul(v, p));
+
+ Float128 r = fputil::quick_add(sum, p);
+
+ return static_cast<double>(r);
+}
+
+} // namespace
+
+LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
+ using FPBits_t = typename fputil::FPBits<double>;
+ constexpr int EXPONENT_BIAS = FPBits_t::EXPONENT_BIAS;
+ constexpr int MANTISSA_WIDTH = FPBits_t::FloatProp::MANTISSA_WIDTH;
+ constexpr uint64_t MANTISSA_MASK = FPBits_t::FloatProp::MANTISSA_MASK;
+ FPBits_t xbits(x);
+ uint64_t x_u = xbits.uintval();
+
+ fputil::DoubleDouble x_dd{0.0, 0.0};
+
+ uint16_t x_exp = xbits.get_unbiased_exponent();
+
+ if (x_exp >= EXPONENT_BIAS) {
+ // |x| >= 1
+ if (LIBC_UNLIKELY(x_u >= 0x4650'0000'0000'0000ULL)) {
+ // x >= 2^102 or x is negative, inf, or NaN
+ if (LIBC_UNLIKELY(x_u > FPBits_t::MAX_NORMAL)) {
+ // x <= -1.0 or x is Inf or NaN
+ if (x_u == 0xbff0'0000'0000'0000ULL) {
+ // x = -1.0
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ return static_cast<double>(FPBits_t::neg_inf());
+ }
+ if (xbits.get_sign() && !xbits.is_nan()) {
+ // x < -1.0
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits_t::build_quiet_nan(0);
+ }
+ // x is +Inf or NaN
+ return x;
+ }
+ x_dd.hi = x;
+ } else {
+ x_dd = fputil::exact_add(x, 1.0);
+ }
+ } else {
+ // |x| < 1
+ if (LIBC_UNLIKELY(xbits.get_unbiased_exponent() <
+ EXPONENT_BIAS - MANTISSA_WIDTH - 1)) {
+ // Quick return when |x| < 2^-53.
+ // Since log(1 + x) = x - x^2/2 + x^3/3 - ...,
+ // for |x| < 2^-53,
+ // x > log(1 + x) > x - x^2 > x(1 - 2^-54) > x - ulp(x)/2
+ // Thus,
+ // log(1 + x) = nextafter(x, -inf) for FE_DOWNWARD, or
+ // FE_TOWARDZERO and x > 0,
+ // = x otherwise.
+ if (LIBC_UNLIKELY(xbits.is_zero()))
+ return x;
+
+ volatile float tp = 1.0f;
+ volatile float tn = -1.0f;
+ bool rdp = (tp - 0x1p-25f != tp);
+ bool rdn = (tn - 0x1p-24f != tn);
+
+ if (x > 0 && rdp) {
+ return FPBits_t(x_u - 1).get_val();
+ }
+
+ if (x < 0 && rdn) {
+ return FPBits_t(x_u + 1).get_val();
+ }
+
+ return x;
+ }
+ x_dd = fputil::exact_add(1.0, x);
+ }
+
+ // At this point, x_dd is the exact sum of 1 + x:
+ // x_dd.hi + x_dd.lo = x + 1.0 exactly.
+ // |x_dd.hi| >= 2^-54
+ // |x_dd.lo| < ulp(x_dd.hi)
+
+ FPBits_t xhi_bits(x_dd.hi);
+ x_u = xhi_bits.uintval();
+ // Range reduction:
+ // Find k such that |x_hi - k * 2^-7| <= 2^-8.
+ int idx = ((x_u & MANTISSA_MASK) + (1ULL << (MANTISSA_WIDTH - 8))) >>
+ (MANTISSA_WIDTH - 7);
+ int x_e = xhi_bits.get_exponent() + (idx >> 7);
+ double e_x = static_cast<double>(x_e);
+
+ // hi is exact
+ // ulp(hi) = ulp(LOG_2_HI) = ulp(LOG_R1_DD[idx].hi) = 2^-43
+ double hi = fputil::multiply_add(e_x, LOG_2_HI, LOG_R1_DD[idx].hi);
+ // lo errors < |e_x| * ulp(LOG_2_LO) + ulp(LOG_R1[idx].lo)
+ // <= 2^11 * 2^(-43-53) = 2^-85
+ double lo = fputil::multiply_add(e_x, LOG_2_LO, LOG_R1_DD[idx].lo);
+
+ // Error bound of e_x * log(2) - log(r1)
+ constexpr double ERR_HI[2] = {0x1.0p-85, 0.0};
+ double err_hi = ERR_HI[hi == 0.0];
+
+ // Scaling factior = 2^(-xh_bits.get_exponent())
+ uint64_t s_u =
+ (static_cast<uint64_t>(EXPONENT_BIAS) << (MANTISSA_WIDTH + 1)) -
+ (x_u & FPBits_t::FloatProp::EXPONENT_MASK);
+ // When the exponent of x is 2^1023, its inverse, 2^(-1023), is subnormal.
+ const double EXPONENT_CORRECTION[2] = {0.0, 0x1.0p-1023};
+ double scaling = FPBits_t(s_u).get_val() + EXPONENT_CORRECTION[s_u == 0];
+ // Normalize arguments:
+ // 1 <= m_dd.hi < 2
+ // |m_dd.lo| < 2^-52.
+ // This is exact.
+ fputil::DoubleDouble m_dd{scaling * x_dd.lo, scaling * x_dd.hi};
+
+ // Perform range reduction:
+ // r * m - 1 = r * (m_dd.hi + m_dd.lo) - 1
+ // = (r * m_dd.hi - 1) + r * m_dd.lo
+ // = v_hi + (v_lo.hi + v_lo.lo)
+ // where:
+ // v_hi = r * m_dd.hi - 1 (exact)
+ // v_lo.hi + v_lo.lo = r * m_dd.lo (exact)
+ // Bounds on the values:
+ // -0x1.69000000000edp-8 < r * m - 1 < 0x1.7f00000000081p-8
+ // |v_lo.hi| <= |r| * |m_dd.lo| < 2^-52
+ // |v_lo.lo| < ulp(v_lo.hi) <= 2^(-52 - 53) = 2^(-105)
+ double r = R1[idx];
+ double v_hi;
+ fputil::DoubleDouble v_lo = fputil::exact_mult(m_dd.lo, r);
+
+ // Perform exact range reduction
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ v_hi = fputil::multiply_add(r, m_dd.hi, -1.0); // Exact.
+#else
+ // c = 1 + idx * 2^-7.
+ double c = FPBits_t((static_cast<uint64_t>(idx) << (MANTISSA_WIDTH - 7)) +
+ uint64_t(0x3FF0'0000'0000'0000ULL))
+ .get_val();
+ v_hi = fputil::multiply_add(r, m_dd.hi - c, RCM1[idx]); // Exact
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+ // Range reduction output:
+ // -0x1.69000000000edp-8 < v_hi + v_lo < 0x1.7f00000000081p-8
+ // |v_dd.lo| < ulp(v_dd.hi) <= 2^(-7 - 53) = 2^-60
+ fputil::DoubleDouble v_dd = fputil::exact_add(v_hi, v_lo.hi);
+ v_dd.lo += v_lo.lo;
+
+ // Exact sum:
+ // r1.hi + r1.lo = e_x * log(2)_hi - log(r)_hi + u
+ fputil::DoubleDouble r1 = fputil::exact_add(hi, v_dd.hi);
+
+ // Overall error is bounded by:
+ // C * ulp(v_sq) + err_hi
+ double v_sq = v_dd.hi * v_dd.hi;
+ double p0 = fputil::multiply_add(v_dd.hi, P_COEFFS[1], P_COEFFS[0]);
+ double p1 = fputil::multiply_add(v_dd.hi, P_COEFFS[3], P_COEFFS[2]);
+ double p2 = fputil::multiply_add(v_dd.hi, P_COEFFS[5], P_COEFFS[4]);
+ double p = fputil::polyeval(v_sq, (v_dd.lo + r1.lo) + lo, p0, p1, p2);
+
+ double err = fputil::multiply_add(v_sq, P_ERR, err_hi);
+
+ double left = r1.hi + (p - err);
+ double right = r1.hi + (p + err);
+
+ // Ziv's test to see if fast pass is accurate enough.
+ if (left == right)
+ return left;
+
+ return log1p_accurate(x_e, idx, v_dd);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/log1p.h b/libc/src/math/log1p.h
new file mode 100644
index 0000000000000..0f25ce6167545
--- /dev/null
+++ b/libc/src/math/log1p.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for log1p -------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_LOG1P_H
+#define LLVM_LIBC_SRC_MATH_LOG1P_H
+
+namespace __llvm_libc {
+
+double log1p(double x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_LOG1P_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 0fff5bde37675..01b8a74463ffd 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1351,6 +1351,20 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+log1p_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ log1p_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.include.math
+ libc.src.math.log1p
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
log1pf_test
NEED_MPFR
diff --git a/libc/test/src/math/log1p_test.cpp b/libc/test/src/math/log1p_test.cpp
new file mode 100644
index 0000000000000..34684329af2d7
--- /dev/null
+++ b/libc/test/src/math/log1p_test.cpp
@@ -0,0 +1,166 @@
+//===-- Unittests for log1p -----------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/log1p.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+using __llvm_libc::testing::tlog;
+
+DECLARE_SPECIAL_CONSTANTS(double)
+
+TEST(LlvmLibcLog1pTest, SpecialNumbers) {
+ EXPECT_FP_EQ(aNaN, __llvm_libc::log1p(aNaN));
+ EXPECT_FP_EQ(inf, __llvm_libc::log1p(inf));
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log1p(neg_inf), FE_INVALID);
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log1p(-2.0), FE_INVALID);
+ EXPECT_FP_EQ(zero, __llvm_libc::log1p(0.0));
+ EXPECT_FP_EQ(neg_zero, __llvm_libc::log1p(-0.0));
+ EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::log1p(-1.0), FE_DIVBYZERO);
+}
+
+TEST(LlvmLibcLog1pTest, TrickyInputs) {
+ constexpr int N = 41;
+ constexpr uint64_t INPUTS[N] = {
+ 0x3ff0000000000000, // x = 1.0
+ 0x4024000000000000, // x = 10.0
+ 0x4059000000000000, // x = 10^2
+ 0x408f400000000000, // x = 10^3
+ 0x40c3880000000000, // x = 10^4
+ 0x40f86a0000000000, // x = 10^5
+ 0x412e848000000000, // x = 10^6
+ 0x416312d000000000, // x = 10^7
+ 0x4197d78400000000, // x = 10^8
+ 0x41cdcd6500000000, // x = 10^9
+ 0x4202a05f20000000, // x = 10^10
+ 0x42374876e8000000, // x = 10^11
+ 0x426d1a94a2000000, // x = 10^12
+ 0x42a2309ce5400000, // x = 10^13
+ 0x42d6bcc41e900000, // x = 10^14
+ 0x430c6bf526340000, // x = 10^15
+ 0x4341c37937e08000, // x = 10^16
+ 0x4376345785d8a000, // x = 10^17
+ 0x43abc16d674ec800, // x = 10^18
+ 0x43e158e460913d00, // x = 10^19
+ 0x4415af1d78b58c40, // x = 10^20
+ 0x444b1ae4d6e2ef50, // x = 10^21
+ 0x4480f0cf064dd592, // x = 10^22
+ 0x3fefffffffef06ad, 0x3fefde0f22c7d0eb, 0x225e7812faadb32f,
+ 0x3fee1076964c2903, 0x3fdfe93fff7fceb0, 0x3ff012631ad8df10,
+ 0x3fefbfdaa448ed98, 0x3fd00a8cefe9a5f8, 0x3fd0b4d870eb22f8,
+ 0x3c90c40cef04efb5, 0x449d2ccad399848e, 0x4aa12ccdffd9d2ec,
+ 0x5656f070b92d36ce, 0x6db06dcb74f76bcc, 0x7f1954e72ffd4596,
+ 0x5671e2f1628093e4, 0x73dac56e2bf1a951, 0x8001bc6879ea14c5,
+ };
+ for (int i = 0; i < N; ++i) {
+ double x = double(FPBits(INPUTS[i]));
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log1p, x,
+ __llvm_libc::log1p(x), 0.5);
+ }
+}
+
+TEST(LlvmLibcLog1pTest, AllExponents) {
+ double x = 0x1.0p-1074;
+ for (int i = -1074; i < 1024; ++i, x *= 2.0) {
+ ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log1p, x,
+ __llvm_libc::log1p(x), 0.5);
+ }
+}
+
+TEST(LlvmLibcLog1pTest, InDoubleRange) {
+ constexpr uint64_t COUNT = 1234561;
+
+ auto test = [&](uint64_t start, uint64_t stop,
+ mpfr::RoundingMode rounding_mode) {
+ mpfr::ForceRoundingMode __r(rounding_mode);
+ uint64_t fails = 0;
+ uint64_t count = 0;
+ uint64_t cc = 0;
+ double mx, mr = 0.0;
+ double tol = 0.5;
+
+ uint64_t step = (stop - start) / COUNT;
+
+ for (uint64_t i = 0, v = start; i <= COUNT; ++i, v += step) {
+ double x = FPBits(v).get_val();
+ if (isnan(x) || isinf(x) || x < 0.0)
+ continue;
+ libc_errno = 0;
+ double result = __llvm_libc::log1p(x);
+ ++cc;
+ if (isnan(result) || isinf(result))
+ continue;
+
+ ++count;
+ // ASSERT_MPFR_MATCH(mpfr::Operation::Log1p, x, result, 0.5);
+ if (!EXPECT_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Log1p, x,
+ result, 0.5, rounding_mode)) {
+ ++fails;
+ while (!EXPECT_MPFR_MATCH_ROUNDING_SILENTLY(
+ mpfr::Operation::Log1p, x, result, tol, rounding_mode)) {
+ mx = x;
+ mr = result;
+ tol *= 2.0;
+ }
+ }
+ }
+ tlog << " Log1p failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
+ if (fails) {
+ EXPECT_MPFR_MATCH(mpfr::Operation::Log1p, mx, mr, 0.5, rounding_mode);
+ }
+ };
+
+ auto test_all_rounding = [&](uint64_t start, uint64_t stop,
+ const char *start_str, const char *stop_str) {
+ tlog << "\n=== Test in range [" << start_str << ", " << stop_str
+ << "] ===\n";
+
+ tlog << "\n Test Rounding To Nearest...\n";
+ test(start, stop, mpfr::RoundingMode::Nearest);
+
+ tlog << "\n Test Rounding Downward...\n";
+ test(start, stop, mpfr::RoundingMode::Downward);
+
+ tlog << "\n Test Rounding Upward...\n";
+ test(start, stop, mpfr::RoundingMode::Upward);
+
+ tlog << "\n Test Rounding Toward Zero...\n";
+ test(start, stop, mpfr::RoundingMode::TowardZero);
+ };
+
+ test_all_rounding(0x0000'0000'0000'0001ULL, 0x0010'0000'0000'0000ULL,
+ "2^-1074", "2^-1022");
+
+ test_all_rounding(0x39B0'0000'0000'0000ULL, 0x3A50'0000'0000'0000ULL,
+ "2^-100", "2^-90");
+
+ test_all_rounding(0x3CD0'0000'0000'0000ULL, 0x3D20'0000'0000'0000ULL, "2^-50",
+ "2^-45");
+
+ test_all_rounding(0x3E10'0000'0000'0000ULL, 0x3E40'0000'0000'0000ULL, "2^-30",
+ "2^-27");
+
+ test_all_rounding(0x3FD0'0000'0000'0000ULL, 0x4010'0000'0000'0000ULL, "0.25",
+ "4.0");
+
+ test_all_rounding(0x4630'0000'0000'0000ULL, 0x4670'0000'0000'0000ULL, "2^100",
+ "2^104");
+
+ test_all_rounding(0x7FD0'0000'0000'0000ULL, 0x7FF0'0000'0000'0000ULL,
+ "2^1022", "2^1024");
+}
More information about the libc-commits
mailing list