diff --git a/libc/config/darwin/arm/entrypoints.txt b/libc/config/darwin/arm/entrypoints.txt --- a/libc/config/darwin/arm/entrypoints.txt +++ b/libc/config/darwin/arm/entrypoints.txt @@ -176,6 +176,7 @@ 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 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -287,6 +287,7 @@ 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 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -292,6 +292,7 @@ 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 --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -169,6 +169,7 @@ 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 --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -408,6 +408,7 @@ FunctionSpec<"log10", RetValSpec, [ArgSpec]>, FunctionSpec<"log10f", RetValSpec, [ArgSpec]>, + FunctionSpec<"log1p", RetValSpec, [ArgSpec]>, FunctionSpec<"log1pf", RetValSpec, [ArgSpec]>, FunctionSpec<"log2", RetValSpec, [ArgSpec]>, diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -114,6 +114,7 @@ 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 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -814,6 +814,26 @@ -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 --- /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(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(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(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(r); +} + +} // namespace + +LLVM_LIBC_FUNCTION(double, log1p, (double x)) { + using FPBits_t = typename fputil::FPBits; + 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(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(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(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(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 --- /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 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -1351,6 +1351,20 @@ 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 --- /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 + +#include +#include + +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(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"); +}