//===-- 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/integer_literals.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY

#include "src/__support/math/common_constants.h"

namespace LIBC_NAMESPACE_DECL {

// 128-bit precision dyadic floating point numbers.
using Float128 = typename fputil::DyadicFloat<128>;

using LIBC_NAMESPACE::operator""_u128;

namespace {

using namespace common_constants_internal;

// 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
alignas(16) 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};

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Extra errors from P is from using x^2 to reduce evaluation latency and
// directional rounding.
constexpr double P_ERR = 0x1.0p-49;

// log(2) with 128-bit precision generated by SageMath with:
// def format_hex(value):
//     l = hex(value)[2:]
//     n = 8
//     x = [l[i:i + n] for i in range(0, len(l), n)]
//     return "0x" + "'".join(x) + "_u128"
// (s, m, e) = RealField(128)(2).log().sign_mantissa_exponent();
// print(format_hex(m));
constexpr Float128 LOG_2(Sign::POS, /*exponent=*/-128, /*mantissa=*/
                         0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128);

// -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("{Sign::POS,", e, ", format_hex(m), "},");
constexpr Float128 LOG_R1[129] = {
    {Sign::POS, 0, 0_u128},
    {Sign::POS, -134, 0x8080abac'46f38946'662d417c'ed007a46_u128},
    {Sign::POS, -133, 0x8102b2c4'9ac23a4f'91d082dc'e3ddcd38_u128},
    {Sign::POS, -133, 0xc2492946'4655f45c'da5f3cc0'b3251dbd_u128},
    {Sign::POS, -132, 0x820aec4f'3a222380'b9e3aea6'c444ef07_u128},
    {Sign::POS, -132, 0xa33576a1'6f1f4c64'521016bd'904dc968_u128},
    {Sign::POS, -132, 0xb3e4a796'a5dac208'27cca0bc'c06c2f92_u128},
    {Sign::POS, -132, 0xd5779687'd887e0d1'a9dda170'56e45ed5_u128},
    {Sign::POS, -132, 0xf7518e00'35c3dd83'606d8909'3278a939_u128},
    {Sign::POS, -131, 0x8cb9de8a'32ab368a'a7c98595'30a45153_u128},
    {Sign::POS, -131, 0x9defad3e'8f73217a'976d3b5b'45f6ca0b_u128},
    {Sign::POS, -131, 0xa6988ae9'03f562ed'3e858f08'597b3a69_u128},
    {Sign::POS, -131, 0xb8069857'560707a3'6a677b4c'8bec22e1_u128},
    {Sign::POS, -131, 0xc99af2ea'ca4c4570'eaf51f66'692844ba_u128},
    {Sign::POS, -131, 0xd273b205'8de1bd49'46bbf837'b4d320c6_u128},
    {Sign::POS, -131, 0xe442c00d'e2591b47'196ab34c'e0bccd12_u128},
    {Sign::POS, -131, 0xed393b1c'22351280'3f4e2e66'0317d55f_u128},
    {Sign::POS, -131, 0xff4489ce'deab2ca6'c17bd40d'8d9291ec_u128},
    {Sign::POS, -130, 0x88bc7411'3f23def1'9c5a0fe3'96f40f1e_u128},
    {Sign::POS, -130, 0x8d515bf1'1fb94f1c'88713268'840cbcc0_u128},
    {Sign::POS, -130, 0x968b0864'3409ceb6'65c0da50'6a088484_u128},
    {Sign::POS, -130, 0x9b2fe580'ac80b17d'411a5b94'4aca8708_u128},
    {Sign::POS, -130, 0xa489ec19'9dab06f2'a9fb6cf0'ecb411b7_u128},
    {Sign::POS, -130, 0xa93f2f25'0dac67d1'cad2fb8d'48054ae0_u128},
    {Sign::POS, -130, 0xadfa035a'a1ed8fdc'149767e4'10316d2c_u128},
    {Sign::POS, -130, 0xb780945b'ab55dce4'34c7bc3d'32750fde_u128},
    {Sign::POS, -130, 0xbc4c6c2a'226399ef'8f6ebcfb'2016a439_u128},
    {Sign::POS, -130, 0xc5f57f59'c7f46155'aa8b6997'a402bf30_u128},
    {Sign::POS, -130, 0xcad2d6e7'b80bf914'2c507fb7'a3d0bf6a_u128},
    {Sign::POS, -130, 0xcfb62038'44b3209a'd0cb02f3'3f79c16c_u128},
    {Sign::POS, -130, 0xd98ec2ba'de71e539'58a98f2a'd65bee9b_u128},
    {Sign::POS, -130, 0xde8439c1'dec56877'4d57da94'5b5d0aaa_u128},
    {Sign::POS, -130, 0xe37fde37'807b84e3'4e9a750b'6b68781d_u128},
    {Sign::POS, -130, 0xe881bf93'2af3dac0'c524848e'3443e040_u128},
    {Sign::POS, -130, 0xf29877ff'38809091'3b020fa1'820c9492_u128},
    {Sign::POS, -130, 0xf7ad6f26'e7ff2ef7'54d2238f'75f969b1_u128},
    {Sign::POS, -130, 0xfcc8e365'9d9bcbec'ca0cdf30'1431b60f_u128},
    {Sign::POS, -129, 0x80f572b1'363487b9'f5bd0b5b'3479d5f4_u128},
    {Sign::POS, -129, 0x86216b3b'0b17188b'163ceae8'8f720f1e_u128},
    {Sign::POS, -129, 0x88bc7411'3f23def1'9c5a0fe3'96f40f1e_u128},
    {Sign::POS, -129, 0x8b5ae65d'67db9acd'f7a51681'26a58b9a_u128},
    {Sign::POS, -129, 0x8dfccb1a'd35ca6ed'5147bdb6'ddcaf59c_u128},
    {Sign::POS, -129, 0x90a22b68'75c6a1f7'ae91aeba'609c8877_u128},
    {Sign::POS, -129, 0x934b1089'a6dc93c1'df5bb3b6'0554e152_u128},
    {Sign::POS, -129, 0x95f783e6'e49a9cfa'4a5004f3'ef063313_u128},
    {Sign::POS, -129, 0x9b5b3bb5'f088b766'd878bbe3'd392be25_u128},
    {Sign::POS, -129, 0x9e1293b9'998c1daa'5b035eae'273a855f_u128},
    {Sign::POS, -129, 0xa0cda11e'af46390d'bb243827'3918db7e_u128},
    {Sign::POS, -129, 0xa38c6e13'8e20d831'f698298a'dddd7f32_u128},
    {Sign::POS, -129, 0xa64f04f0'b961df76'e4f5275c'2d15c21f_u128},
    {Sign::POS, -129, 0xa9157039'c51ebe70'8164c759'686a2209_u128},
    {Sign::POS, -129, 0xabdfba9e'468fd6f6'f72ea077'49ce6bd3_u128},
    {Sign::POS, -129, 0xaeadeefa'caf97d35'7dd6e688'ebb13b03_u128},
    {Sign::POS, -129, 0xb1801859'd56249dc'18ce51ff'f99479cd_u128},
    {Sign::POS, -129, 0xb45641f4'e350a0d3'2756eba0'0bc33978_u128},
    {Sign::POS, -129, 0xb7307735'78cb90b2'be1116c3'466beb6d_u128},
    {Sign::POS, -129, 0xba0ec3b6'33dd8b09'49dc60b2'b059a60b_u128},
    {Sign::POS, -129, 0xbcf13343'e7d9ec7d'2efd1778'1bb3afec_u128},
    {Sign::POS, -129, 0xbfd7d1de'c0a8df6f'37eda996'244bccb0_u128},
    {Sign::POS, -129, 0xc2c2abbb'6e5fd56f'33337789'd592e296_u128},
    {Sign::POS, -129, 0xc5b1cd44'596fa51e'1a18fb8f'9f9ef280_u128},
    {Sign::POS, -129, 0xc8a5431a'dfb44ca5'688ce7c1'a75e341a_u128},
    {Sign::POS, -129, 0xcb9d1a18'9ab56e76'2d7e9307'c70c0668_u128},
    {Sign::POS, -129, 0xcb9d1a18'9ab56e76'2d7e9307'c70c0668_u128},
    {Sign::POS, -129, 0xce995f50'af69d861'ef2f3f4f'861ad6a9_u128},
    {Sign::POS, -129, 0xd19a2011'27d3c645'7f9d79f5'1dcc7301_u128},
    {Sign::POS, -129, 0xd49f69e4'56cf1b79'5f53bd2e'406e66e7_u128},
    {Sign::POS, -129, 0xd7a94a92'466e833a'ad88bba7'd0cee8e0_u128},
    {Sign::POS, -129, 0xdab7d022'31484a92'96c20cca'6efe2ac5_u128},
    {Sign::POS, -129, 0xddcb08dc'0717d85b'f40a666c'87842843_u128},
    {Sign::POS, -129, 0xe0e30349'fd1cec80'7fe8e180'2aba24d6_u128},
    {Sign::POS, -129, 0xe0e30349'fd1cec80'7fe8e180'2aba24d6_u128},
    {Sign::POS, -129, 0xe3ffce3a'2aa64922'3eadb651'b49ac53a_u128},
    {Sign::POS, -129, 0xe72178c0'323a1a0f'304e1653'e71d9973_u128},
    {Sign::POS, -129, 0xea481236'f7d35baf'e9a767a8'0d6d97e8_u128},
    {Sign::POS, -129, 0xed73aa42'64b0ade9'4f91cf4b'33e42998_u128},
    {Sign::POS, -129, 0xed73aa42'64b0ade9'4f91cf4b'33e42998_u128},
    {Sign::POS, -129, 0xf0a450d1'39366ca6'fc66eb64'08ff6433_u128},
    {Sign::POS, -129, 0xf3da161e'ed6b9aaf'ac8d42f7'8d3e65d3_u128},
    {Sign::POS, -129, 0xf7150ab5'a09f27f4'5a470250'd40ebe90_u128},
    {Sign::POS, -129, 0xf7150ab5'a09f27f4'5a470250'd40ebe90_u128},
    {Sign::POS, -129, 0xfa553f70'18c966f2'b780a545'a1b54dcf_u128},
    {Sign::POS, -129, 0xfd9ac57b'd244217e'8f05924d'258c14c5_u128},
    {Sign::POS, -128, 0x8072d72d'903d588b'89d1b09c'70c4010a_u128},
    {Sign::POS, -128, 0x8072d72d'903d588b'89d1b09c'70c4010a_u128},
    {Sign::POS, -128, 0x821b05f3'b01d6774'030d58c3'f7e2ea1f_u128},
    {Sign::POS, -128, 0x83c5f829'9e2b4091'20f6fafe'8fbb68b9_u128},
    {Sign::POS, -128, 0x8573b716'82a7d21a'e21f9f89'c1ab80b2_u128},
    {Sign::POS, -128, 0x8573b716'82a7d21a'e21f9f89'c1ab80b2_u128},
    {Sign::POS, -128, 0x87244c30'8e670a66'01e005d0'6dbfa8f8_u128},
    {Sign::POS, -128, 0x88d7c11e'3ad53cdc'223111a7'07b6de2c_u128},
    {Sign::POS, -128, 0x88d7c11e'3ad53cdc'223111a7'07b6de2c_u128},
    {Sign::POS, -128, 0x8a8e1fb7'94b09134'2eb628db'a173c82d_u128},
    {Sign::POS, -128, 0x8c477207'91e53313'be2ad194'15fe25a5_u128},
    {Sign::POS, -128, 0x8c477207'91e53313'be2ad194'15fe25a5_u128},
    {Sign::POS, -128, 0x8e03c24d'73003959'bddae1cc'ce247838_u128},
    {Sign::POS, -128, 0x8fc31afe'30b2c6de'9b00bf16'7e95da67_u128},
    {Sign::POS, -128, 0x8fc31afe'30b2c6de'9b00bf16'7e95da67_u128},
    {Sign::POS, -128, 0x918586c5'f5e4bf01'9b92199e'd1a4bab1_u128},
    {Sign::POS, -128, 0x934b1089'a6dc93c1'df5bb3b6'0554e152_u128},
    {Sign::POS, -128, 0x934b1089'a6dc93c1'df5bb3b6'0554e152_u128},
    {Sign::POS, -128, 0x9513c368'76083695'f3cbc416'a2418012_u128},
    {Sign::POS, -128, 0x96dfaabd'86fa1646'be1188fb'c94e2f15_u128},
    {Sign::POS, -128, 0x96dfaabd'86fa1646'be1188fb'c94e2f15_u128},
    {Sign::POS, -128, 0x98aed221'a03458b6'1d2f8932'1647b358_u128},
    {Sign::POS, -128, 0x98aed221'a03458b6'1d2f8932'1647b358_u128},
    {Sign::POS, -128, 0x9a81456c'ec642e0f'e549f9aa'ea3cb5e1_u128},
    {Sign::POS, -128, 0x9c5710b8'cbb73a42'a2554b2d'd4619e63_u128},
    {Sign::POS, -128, 0x9c5710b8'cbb73a42'a2554b2d'd4619e63_u128},
    {Sign::POS, -128, 0x9e304061'b5fda919'30603d87'b6df81ad_u128},
    {Sign::POS, -128, 0x9e304061'b5fda919'30603d87'b6df81ad_u128},
    {Sign::POS, -128, 0xa00ce109'2e5498c3'67879c5a'30cd1242_u128},
    {Sign::POS, -128, 0xa00ce109'2e5498c3'67879c5a'30cd1242_u128},
    {Sign::POS, -128, 0xa1ecff97'c91e267b'0b7efae0'8e597e16_u128},
    {Sign::POS, -128, 0xa3d0a93f'45169a4a'83594fab'088c0d65_u128},
    {Sign::POS, -128, 0xa3d0a93f'45169a4a'83594fab'088c0d65_u128},
    {Sign::POS, -128, 0xa5b7eb7c'b860fb88'af6a62a0'dec6e073_u128},
    {Sign::POS, -128, 0xa5b7eb7c'b860fb88'af6a62a0'dec6e073_u128},
    {Sign::POS, -128, 0xa7a2d41a'd270c9d7'49362382'a768847a_u128},
    {Sign::POS, -128, 0xa7a2d41a'd270c9d7'49362382'a768847a_u128},
    {Sign::POS, -128, 0xa9917134'33c2b998'8ba4aea6'14d05701_u128},
    {Sign::POS, -128, 0xa9917134'33c2b998'8ba4aea6'14d05701_u128},
    {Sign::POS, -128, 0xab83d135'dc633301'7fe6607b'a902ef3c_u128},
    {Sign::POS, -128, 0xab83d135'dc633301'7fe6607b'a902ef3c_u128},
    {Sign::POS, -128, 0xad7a02e1'b24efd31'd60864fd'949b4bd3_u128},
    {Sign::POS, -128, 0xad7a02e1'b24efd31'd60864fd'949b4bd3_u128},
    {Sign::POS, -128, 0xaf741551'20c9011c'066d235e'e63073dd_u128},
    {Sign::POS, -128, 0xaf741551'20c9011c'066d235e'e63073dd_u128},
    {Sign::POS, 0, 0_u128},
};

// 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("{Sign::POS," if (s == -1) else "{Sign::NEG,", e, ",
//         format_hex(m), "},");
constexpr Float128 LOG_R2[198] = {
    {Sign::NEG, -135, 0xb67dab2a'1a5742a4'a0e061c5'f7431c5e_u128},
    {Sign::NEG, -135, 0xb4807f24'af682939'5d5bfe7b'969ed6ec_u128},
    {Sign::NEG, -135, 0xb2834b35'b4d54d5f'4d08702d'dfabc23f_u128},
    {Sign::NEG, -135, 0xb0860f5c'eba9be95'd4d36650'8b9953df_u128},
    {Sign::NEG, -135, 0xae68f71a'a09e8847'ac18a289'f8f214a9_u128},
    {Sign::NEG, -135, 0xac6baaee'd676e8f1'd5b42054'abb88c45_u128},
    {Sign::NEG, -135, 0xaa6e56d8'7cd632d6'09809d58'ee484964_u128},
    {Sign::NEG, -135, 0xa870fad7'54bb8791'b9e6fc7c'72f06d73_u128},
    {Sign::NEG, -135, 0xa67396eb'1f231892'6f78d6d0'105c00e2_u128},
    {Sign::NEG, -135, 0xa4762b13'9d0626e7'028f7126'29209148_u128},
    {Sign::NEG, -135, 0xa258dfd1'0aedaa67'c98d898e'f172df02_u128},
    {Sign::NEG, -135, 0xa05b63a3'73e60a83'fcc37c3c'3062bfa1_u128},
    {Sign::NEG, -135, 0x9e5ddf89'cf42f501'3eb450db'05763c36_u128},
    {Sign::NEG, -135, 0x9c605383'ddf1b88c'7146a86f'd458b775_u128},
    {Sign::NEG, -135, 0x9a62bf91'60dcb286'c20a0c92'81474436_u128},
    {Sign::NEG, -135, 0x986523b2'18eb4ed6'cdc57316'ec4aebc3_u128},
    {Sign::NEG, -135, 0x96677fe5'c70207b9'c060dad7'4cef4273_u128},
    {Sign::NEG, -135, 0x9449f92d'2ff44633'ed8def1a'3e433499_u128},
    {Sign::NEG, -135, 0x924c4507'3220b5e0'3ce7a1f8'5c27b4fc_u128},
    {Sign::NEG, -135, 0x904e88f3'68fea63f'f2ca8934'49f7f2cb_u128},
    {Sign::NEG, -135, 0x8e50c4f1'956699ed'8d77d9fa'bd2853cf_u128},
    {Sign::NEG, -135, 0x8c52f901'782e20ec'93e828d7'5b58ded4_u128},
    {Sign::NEG, -135, 0x8a552522'd227d87a'9f9605b0'53c5acf0_u128},
    {Sign::NEG, -135, 0x88574955'64236ae0'62a14939'3bca7241_u128},
    {Sign::NEG, -135, 0x86398719'b66bac7c'aea6b56c'e89203d4_u128},
    {Sign::NEG, -135, 0x843b9aef'044e4dcc'0242bd86'd00609b2_u128},
    {Sign::NEG, -135, 0x823da6d4'c89c6927'daabf927'74bac84e_u128},
    {Sign::NEG, -135, 0x803faaca'c419abf2'a1c6f3fc'242ef8d0_u128},
    {Sign::NEG, -136, 0xfc834da1'6f0d9f57'a225ebc0'2e6d9dd4_u128},
    {Sign::NEG, -136, 0xf88735cc'c7433381'c33f6ad3'40ae18a9_u128},
    {Sign::NEG, -136, 0xf48b0e17'1249b6bc'70b2a4d3'8a242244_u128},
    {Sign::NEG, -136, 0xf08ed67f'd190e280'1d548190'48b811b0_u128},
    {Sign::NEG, -136, 0xec52ca07'ed95f236'9c21b650'afe9ede0_u128},
    {Sign::NEG, -136, 0xe85671ad'ecd28aac'935519c9'6d30e463_u128},
    {Sign::NEG, -136, 0xe45a0970'dc912ca7'ba88f6f2'e2672cfe_u128},
    {Sign::NEG, -136, 0xe05d9150'3e298bc8'0b1a8b84'657ae069_u128},
    {Sign::NEG, -136, 0xdc61094b'92ed70ef'ea3bff8d'197b20a1_u128},
    {Sign::NEG, -136, 0xd8647162'5c28b9e5'cdbb931d'6fecc249_u128},
    {Sign::NEG, -136, 0xd467c994'1b2158f5'd971d560'd5f00820_u128},
    {Sign::NEG, -136, 0xd06b11e0'51175493'75563561'244c090b_u128},
    {Sign::NEG, -136, 0xcc6e4a46'7f44c6fa'dc393c9a'3f3b380f_u128},
    {Sign::NEG, -136, 0xc831a4c6'f6fa709d'e6abe6e9'e4ee2096_u128},
    {Sign::NEG, -136, 0xc434bc61'24a0f16e'3ce3c822'8583a66e_u128},
    {Sign::NEG, -136, 0xc037c413'c61bfd93'b96a79f5'c5a4963a_u128},
    {Sign::NEG, -136, 0xbc3abbde'5c8d9bde'aaef2733'7008679f_u128},
    {Sign::NEG, -136, 0xb83da3c0'6911e509'a49a3fca'ddc8bc5a_u128},
    {Sign::NEG, -136, 0xb4407bb9'6cbf035a'e0254feb'785362fa_u128},
    {Sign::NEG, -136, 0xb04343c8'e8a53245'9893a4e2'5ab9dc95_u128},
    {Sign::NEG, -136, 0xac45fbee'5dcebe0b'5d8b0f40'a3708915_u128},
    {Sign::NEG, -136, 0xa848a429'4d40035d'5f4c11c2'c7a58c69_u128},
    {Sign::NEG, -136, 0xa44b3c79'37f76efd'b348cc5d'f706ffba_u128},
    {Sign::NEG, -136, 0xa04dc4dd'9eed7d60'9159f2c5'5a18befd_u128},
    {Sign::NEG, -136, 0x9c106456'3058bef3'bdfdee41'fe6a5a02_u128},
    {Sign::NEG, -136, 0x9812cbe3'46475a24'4580ddf8'9853254d_u128},
    {Sign::NEG, -136, 0x94152383'53489ffb'ac75e10d'61fc3ee8_u128},
    {Sign::NEG, -136, 0x90176b35'd83ce8e2'cad9b30b'29736155_u128},
    {Sign::NEG, -136, 0x8c19a2fa'55fe9b14'6f881deb'98fc45f3_u128},
    {Sign::NEG, -136, 0x881bcad0'4d622a3e'70a04b63'b7248c96_u128},
    {Sign::NEG, -136, 0x841de2b7'3f361722'b4823fb4'8035eddd_u128},
    {Sign::NEG, -136, 0x801feaae'ac42ef38'3364ccb5'b13cd47f_u128},
    {Sign::NEG, -137, 0xf843c56c'2a969897'e306977b'049f0ad5_u128},
    {Sign::NEG, -137, 0xf0479599'f617a843'e3c4d9e9'619bc045_u128},
    {Sign::NEG, -137, 0xe84b45e5'bc76702c'4356d525'b5e6432d_u128},
    {Sign::NEG, -137, 0xe04ed64e'7f14697a'7839dcd7'989339ab_u128},
    {Sign::NEG, -137, 0xd85246d3'3f47230b'4e21f045'ecb76f23_u128},
    {Sign::NEG, -137, 0xd0559772'fe5840b0'902e248d'd4ba9b28_u128},
    {Sign::NEG, -137, 0xc858c82c'bd857a72'a4444906'7ef92e01_u128},
    {Sign::NEG, -137, 0xc05bd8ff'7e009bd2'17926207'cc22e4e6_u128},
    {Sign::NEG, -137, 0xb85ec9ea'40ef8309'1c349622'f3fa5d82_u128},
    {Sign::NEG, -137, 0xafe1c6ec'e1a058dd'97fa2fd0'c9dc723e_u128},
    {Sign::NEG, -137, 0xa7e47606'048b1a65'983e8089'7cf1e60f_u128},
    {Sign::NEG, -137, 0x9fe70534'1d236102'7199cd06'ae5d39b3_u128},
    {Sign::NEG, -137, 0x97e97476'2c5e8f58'43cd18a7'2a051a96_u128},
    {Sign::NEG, -137, 0x8febc3cb'332616ff'7b6d1248'c3e1fd40_u128},
    {Sign::NEG, -137, 0x87edf332'325777c5'f5572a88'14c703af_u128},
    {Sign::NEG, -138, 0xffe00554'55887de0'26828c92'649a3a39_u128},
    {Sign::NEG, -138, 0xefe3e464'3a640cf3'82c550bd'1216d82a_u128},
    {Sign::NEG, -138, 0xdfe78392'14b4e8ae'da6959f7'f0e01bf0_u128},
    {Sign::NEG, -138, 0xcfeae2db'e5d6736d'da93e2fa'85a8f214_u128},
    {Sign::NEG, -138, 0xbfee023f'af0c2480'b47505bf'a5a03b06_u128},
    {Sign::NEG, -138, 0xaff0e1bb'718186ad'b1475a51'80a43520_u128},
    {Sign::NEG, -138, 0x9ff3814d'2e4a36b2'a8740b91'c95df537_u128},
    {Sign::NEG, -138, 0x8ff5e0f2'e661e1c6'57d895d3'5921b59c_u128},
    {Sign::NEG, -139, 0xfff00155'35588833'3c56c598'c659c2a3_u128},
    {Sign::NEG, -139, 0xdff3c0e4'97ea4eb1'2ef8ec33'ed9d782a_u128},
    {Sign::NEG, -139, 0xbff7008f'f5e0c257'379eba7e'6465ff63_u128},
    {Sign::NEG, -139, 0x9ff9c053'5073a370'3f972b78'3fcab757_u128},
    {Sign::NEG, -140, 0xfff80055'51558885'de026e27'1ee0549d_u128},
    {Sign::NEG, -140, 0xbffb8023'febc0c25'eceb47ea'01f6c632_u128},
    {Sign::NEG, -141, 0xfffc0015'54d55888'7333c578'57e1ed52_u128},
    {Sign::NEG, -142, 0xfffe0005'55455588'87dde026'fa704374_u128},
    {Sign::POS, 0, 0_u128},
    {Sign::POS, -141, 0x80010002'aab2aac4'44999abe'2fe2cc65_u128},
    {Sign::POS, -140, 0x8002000a'aaeaac44'4eef3815'81464ccb_u128},
    {Sign::POS, -140, 0xc0048024'01440c26'dfeb4850'85f6f454_u128},
    {Sign::POS, -139, 0x8004002a'acaac445'99abe3be'3a1c6e93_u128},
    {Sign::POS, -139, 0xa0064053'5a37a37a'6bc1e20e'ac8448b4_u128},
    {Sign::POS, -139, 0xc0090090'0a20c275'979eedc0'64c242fd_u128},
    {Sign::POS, -139, 0xe00c40e4'bd6e4efd'c72446cc'1bf728bd_u128},
    {Sign::POS, -138, 0x800800aa'baac446e'f381b821'bbb569e5_u128},
    {Sign::POS, -138, 0x900a20f3'19a3e273'569b26aa'a485ea5c_u128},
    {Sign::POS, -138, 0xa00c814d'7c6a37f8'2dcf56c8'3c80b028_u128},
    {Sign::POS, -138, 0xb00f21bb'e3e388ee'5f697682'84463b9b_u128},
    {Sign::POS, -138, 0xc0120240'510c284c'b48ea6c0'5e2773a1_u128},
    {Sign::POS, -138, 0xd01522dc'c4f87991'14d9d761'96d8043a_u128},
    {Sign::POS, -138, 0xe0188393'40d4f241'e016a611'a4415d72_u128},
    {Sign::POS, -138, 0xf01c2465'c5e61b6f'661e135f'49a47c40_u128},
    {Sign::POS, -137, 0x801002ab'2ac4499a'be6bf0fa'435e8383_u128},
    {Sign::POS, -137, 0x88121333'7898871e'9a31ba0c'bc030353_u128},
    {Sign::POS, -137, 0x901443cc'cd362c9f'54b57dfe'0c4c840f_u128},
    {Sign::POS, -137, 0x98169478'296fad41'7ad1e9c3'15328f7e_u128},
    {Sign::POS, -137, 0xa0190536'8e2389b3'1f3f686c'f3d6be22_u128},
    {Sign::POS, -137, 0xa81b9608'fc3c50ec'f105b66e'c4703ede_u128},
    {Sign::POS, -137, 0xb01e46f0'74b0a0f3'610848c6'8df4d233_u128},
    {Sign::POS, -137, 0xb7a0e9ed'7613acb0'2e0efddf'33a20464_u128},
    {Sign::POS, -137, 0xbfa3d900'8e042ffb'c2cdb3c7'50f127b4_u128},
    {Sign::POS, -137, 0xc7a6e82b'a36a7073'bd953378'6d3f4c49_u128},
    {Sign::POS, -137, 0xcfaa176f'b76c8eb1'82e237c9'a4d450e3_u128},
    {Sign::POS, -137, 0xd7ad66cd'cb3cbe14'c00b46a4'd0e3dfd0_u128},
    {Sign::POS, -137, 0xdfb0d646'e0194584'ea999c0d'f8546710_u128},
    {Sign::POS, -137, 0xe7b465db'f74c8032'cec6c2a9'ad974f4f_u128},
    {Sign::POS, -137, 0xefb8158e'122cde5a'2d2045da'1570a07c_u128},
    {Sign::POS, -137, 0xf7bbe55e'321ce603'6752e9b2'381e3edc_u128},
    {Sign::POS, -137, 0xffbfd54d'588b33c5'3c1ed527'28e00e40_u128},
    {Sign::POS, -136, 0x83e1f2ae'43793dc3'493b0d87'3fb9a340_u128},
    {Sign::POS, -136, 0x87e40ac6'5f6cc4a0'29e38750'c9d26893_u128},
    {Sign::POS, -136, 0x8be632ef'80e9a0df'aab9e832'7258ac3f_u128},
    {Sign::POS, -136, 0x8fe86b2a'28bf51b3'28bc403d'8a5f3c63_u128},
    {Sign::POS, -136, 0x93eab376'd7c36377'f720c1c9'7227fcdc_u128},
    {Sign::POS, -136, 0x97ed0bd6'0ed17018'6ad9a3e3'd11b66c1_u128},
    {Sign::POS, -136, 0x9bef7448'4ecb1f6c'edb27b79'c90b4019_u128},
    {Sign::POS, -136, 0x9fb1c4cd'27012e19'a092a0d7'ab21722a_u128},
    {Sign::POS, -136, 0xa3b44c65'b71c2d85'535d52f0'939a4d02_u128},
    {Sign::POS, -136, 0xa7b6e412'cadcb3dc'90a57e11'edc1864e_u128},
    {Sign::POS, -136, 0xabb98bd4'e33c4381'68e9c901'60031159_u128},
    {Sign::POS, -136, 0xafbc43ac'813a6ea3'bf60594f'929adeb8_u128},
    {Sign::POS, -136, 0xb3bf0b9a'25dcd7a2'8a421588'86775205_u128},
    {Sign::POS, -136, 0xb7c1e39e'522f316d'1ab45417'663dee9e_u128},
    {Sign::POS, -136, 0xbbc4cbb9'87433fe4'6c51ae3c'e1aea68a_u128},
    {Sign::POS, -136, 0xbfc7c3ec'4630d83c'7c52ae8b'40ebabb7_u128},
    {Sign::POS, -136, 0xc3cacc37'1015e15d'a857126f'7cfaaa67_u128},
    {Sign::POS, -136, 0xc7cde49a'66165446'14d05662'cd29464a_u128},
    {Sign::POS, -136, 0xcb90da16'44d29bb7'8379db06'ef3cd6bb_u128},
    {Sign::POS, -136, 0xcf9411aa'99ddb7de'9025f4c6'7dd38bb6_u128},
    {Sign::POS, -136, 0xd3975958'f681086d'd6f8a61c'892032ee_u128},
    {Sign::POS, -136, 0xd79ab121'dbf8714c'9a2f20b4'e2332d47_u128},
    {Sign::POS, -136, 0xdb9e1905'cb85ea59'3c767d61'f51d375b_u128},
    {Sign::POS, -136, 0xdfa19105'46717fca'd4b2bd65'bb25493c_u128},
    {Sign::POS, -136, 0xe3a51920'ce095292'c96c1254'a30ef91f_u128},
    {Sign::POS, -136, 0xe7a8b158'e3a198be'73e324ce'0946b214_u128},
    {Sign::POS, -136, 0xebac59ae'08949dd8'cacd125a'12bac62c_u128},
    {Sign::POS, -136, 0xef6fd620'b2b7a503'cafdc272'27b71eaa_u128},
    {Sign::POS, -136, 0xf3739daf'959aaafc'688d4282'f6026aa3_u128},
    {Sign::POS, -136, 0xf777755d'03f4e0b6'e54e9e38'04464cdd_u128},
    {Sign::POS, -136, 0xfb7b5d29'7f388a12'cb78b383'f4b59dce_u128},
    {Sign::POS, -136, 0xff7f5515'88de024f'ee055fc5'15062c04_u128},
    {Sign::POS, -135, 0x81c1ae90'd131de38'207812b4'3382acdd_u128},
    {Sign::POS, -135, 0x83c3baa7'26a721cc'dc90c4c4'b61f3a87_u128},
    {Sign::POS, -135, 0x85c5cece'05941dbc'1a03f13f'b2c978b1_u128},
    {Sign::POS, -135, 0x87c7eb05'aec1304f'b36f282e'83a7dc36_u128},
    {Sign::POS, -135, 0x89a9eccd'56a980c0'd82a4661'6d4c393f_u128},
    {Sign::POS, -135, 0x8bac18a6'40185360'bc6ff847'13c9babd_u128},
    {Sign::POS, -135, 0x8dae4c90'b22574f4'9f7942a5'16fc2d8a_u128},
    {Sign::POS, -135, 0x8fb0888c'eda546ab'15e50cfd'9b29b427_u128},
    {Sign::POS, -135, 0x91b2cc9b'336f3718'9f465296'ae7dd49a_u128},
    {Sign::POS, -135, 0x93b518bb'c45dc268'b49c1eb9'b348e6e4_u128},
    {Sign::POS, -135, 0x95b76cee'e14e728e'daa320cd'64c9d9c7_u128},
    {Sign::POS, -135, 0x9799a333'de49b963'75a91950'ffe1e3b5_u128},
    {Sign::POS, -135, 0x999c070b'a32068cd'5c6abcbf'43f03f14_u128},
    {Sign::POS, -135, 0x9b9e72f6'b295ad4f'5a9e7f26'5d1ed157_u128},
    {Sign::POS, -135, 0x9da0e6f5'4d9318fd'efeb98d0'2a195c17_u128},
    {Sign::POS, -135, 0x9fa36307'b5054ca8'2aa503a3'110ab5a7_u128},
    {Sign::POS, -135, 0xa1a5e72e'29dbf808'd0fe7e05'869eb825_u128},
    {Sign::POS, -135, 0xa3884a68'a750cb10'e80a28f4'e1e500d2_u128},
    {Sign::POS, -135, 0xa58ade36'aeef9f0b'53106415'1ca6e30b_u128},
    {Sign::POS, -135, 0xa78d7a19'82c4b08f'27c01ffa'8e2e3c4b_u128},
    {Sign::POS, -135, 0xa9901e11'63cbbbf5'7ba9408d'c857d568_u128},
    {Sign::POS, -135, 0xab92ca1e'93038d76'104d1e33'31d3b4fa_u128},
    {Sign::POS, -135, 0xad957e41'516e0158'9343c846'fcdf9137_u128},
    {Sign::POS, -135, 0xaf780e79'b2514889'3977e89a'ec59bfa2_u128},
    {Sign::POS, -135, 0xb17ad246'ef3713bc'913d4e3d'c55c3e6e_u128},
    {Sign::POS, -135, 0xb37d9e2a'7a56b09d'777b52a9'e70d8bcc_u128},
    {Sign::POS, -135, 0xb5807224'94be0c91'55de916f'd30591de_u128},
    {Sign::POS, -135, 0xb7834e35'7f7e2600'e79cfb37'be2861e4_u128},
    {Sign::POS, -135, 0xb986325d'7bab0c89'90983104'd3805389_u128},
    {Sign::POS, -135, 0xbb68ef9c'254aa378'59e3b2ec'71ce64f4_u128},
    {Sign::POS, -135, 0xbd6be371'8c77636f'e83183bf'3dd612ef_u128},
    {Sign::POS, -135, 0xbf6edf5e'c44d9d35'c4e3b0ac'2fd52b7f_u128},
};

// 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("{Sign::POS," if (s == -1) else "{Sign::NEG,", e, ",
//         format_hex(m), "},");
constexpr Float128 LOG_R3[139] = {
    {Sign::NEG, -142, 0x89ff6b38'd5de2622'e39d3faf'42340ed7_u128},
    {Sign::NEG, -142, 0x87ff6f80'ccb40f16'7ff33266'82c02485_u128},
    {Sign::NEG, -142, 0x85ff73b8'c3cdf731'5caf4fbe'343cf928_u128},
    {Sign::NEG, -142, 0x83ff77e0'bb2ade79'cdb6e554'348f7fe8_u128},
    {Sign::NEG, -142, 0x81ff7bf8'b2c9c4f6'0ef009c2'457de25d_u128},
    {Sign::NEG, -143, 0xffff0001'55535558'8883333c'57b57c74_u128},
    {Sign::NEG, -143, 0xfbff07f1'45931f44'f32668f3'9c70d183_u128},
    {Sign::NEG, -143, 0xf7ff0fc1'3650e7bd'459a73c6'a6486fe3_u128},
    {Sign::NEG, -143, 0xf3ff1771'278aaecd'37b18cca'7dd3a29f_u128},
    {Sign::NEG, -143, 0xefff1f01'193e7480'513f610d'21bcfc78_u128},
    {Sign::NEG, -143, 0xebff2671'0b6a38e1'ea190b95'c0690b7b_u128},
    {Sign::NEG, -143, 0xe7ff2dc0'fe0bfbfd'2a150f64'f0ad1743_u128},
    {Sign::NEG, -143, 0xe3ff34f0'f121bddd'090b5174'e995e9d1_u128},
    {Sign::NEG, -143, 0xdfff3c00'e4a97e8c'4ed512b9'b93ea2bf_u128},
    {Sign::NEG, -143, 0xdbff42f0'd8a13e15'934cea21'7ab794a2_u128},
    {Sign::NEG, -143, 0xd7ff49c0'cd06fc83'3e4ebe94'8afd2c76_u128},
    {Sign::NEG, -143, 0xd3ff5070'c1d8b9df'87b7c0f5'bcfee2e1_u128},
    {Sign::NEG, -143, 0xcfff5700'b7147634'77666622'8cb6371b_u128},
    {Sign::NEG, -143, 0xcbff5d70'acb8318b'e53a60f3'514db358_u128},
    {Sign::NEG, -143, 0xc7ff63c0'a2c1ebef'79149c3b'6e57fa86_u128},
    {Sign::NEG, -143, 0xc3ff69f0'992fa568'aad734c9'8416df2a_u128},
    {Sign::NEG, -143, 0xbfff7000'8fff5e00'c2657367'9ed28334_u128},
    {Sign::NEG, -143, 0xbbff75f0'872f15c0'd7a3c6db'6540809f_u128},
    {Sign::NEG, -143, 0xb7ff7bc0'7ebcccb1'd277bde6'45fb1aad_u128},
    {Sign::NEG, -143, 0xb3ff8170'76a682dc'6ac80145'a4087793_u128},
    {Sign::NEG, -143, 0xafff8700'6eea3849'287c4db3'0271e265_u128},
    {Sign::NEG, -143, 0xabff8c70'6785ed00'637d6de4'2eeb151e_u128},
    {Sign::NEG, -143, 0xa7ff91c0'6077a10a'43b5348b'6b898a8c_u128},
    {Sign::NEG, -143, 0xa3ff96f0'59bd546e'c10e7657'978bd7f6_u128},
    {Sign::NEG, -143, 0x9fff9c00'53550735'a37503f4'57310e59_u128},
    {Sign::NEG, -143, 0x9bffa0f0'4d3cb966'82d5a40a'3aa022ff_u128},
    {Sign::NEG, -143, 0x97ffa5c0'47726b08'c71e0d3e'e3df5f4d_u128},
    {Sign::NEG, -143, 0x93ffaa70'41f41c23'a83ce035'2bdbd79b_u128},
    {Sign::NEG, -143, 0x8fffaf00'3cbfccbe'2e21a18d'4680e8e4_u128},
    {Sign::NEG, -143, 0x8bffb370'37d37cdf'30bcb3e4'e5dfbd28_u128},
    {Sign::NEG, -143, 0x87ffb7c0'332d2c8d'57ff51d7'5c66d64a_u128},
    {Sign::NEG, -143, 0x83ffbbf0'2ecadbcf'1bdb87fd'be299f43_u128},
    {Sign::NEG, -144, 0xffff8000'55551555'88885dde'02700703_u128},
    {Sign::NEG, -144, 0xf7ff87e0'4d94724c'd259ca80'3a0c1870_u128},
    {Sign::NEG, -144, 0xefff8f80'464fce8f'e5141308'51c7070a_u128},
    {Sign::NEG, -144, 0xe7ff96e0'3f832a2a'30a16898'f3073a64_u128},
    {Sign::NEG, -144, 0xdfff9e00'392a8526'c4ed6451'7b2949ce_u128},
    {Sign::NEG, -144, 0xd7ffa4e0'3341df90'51e4fb4e'32cf6350_u128},
    {Sign::NEG, -144, 0xcfffab80'2dc53971'277672a8'8350bcce_u128},
    {Sign::NEG, -144, 0xc7ffb1e0'28b092d3'35915377'2a490f06_u128},
    {Sign::NEG, -144, 0xbfffb800'23ffebc0'0c265ece'6b481a0e_u128},
    {Sign::NEG, -144, 0xb7ffbde0'1faf4440'db2781c0'3fa132f6_u128},
    {Sign::NEG, -144, 0xafffc380'1bba9c5e'7287c95c'845ada33_u128},
    {Sign::NEG, -144, 0xa7ffc8e0'181df421'423b56b1'263e5a77_u128},
    {Sign::NEG, -144, 0x9fffce00'14d54b91'5a3752ca'4c076fa3_u128},
    {Sign::NEG, -144, 0x97ffd2e0'11dca2b6'6a71e2b2'7eb3f573_u128},
    {Sign::NEG, -144, 0x8fffd780'0f2ff997'c2e21b72'cff39d8f_u128},
    {Sign::NEG, -144, 0x87ffdbe0'0ccb503c'537ff612'feb7ac9e_u128},
    {Sign::NEG, -145, 0xffffc000'15554d55'58888733'33c57c18_u128},
    {Sign::NEG, -145, 0xefffc7c0'1193f9d1'fa514218'42311c42_u128},
    {Sign::NEG, -145, 0xdfffcf00'0e4aa5fa'2c4ed6de'475b942c_u128},
    {Sign::NEG, -145, 0xcfffd5c0'0b7151d8'ce77678c'bb6fcb88_u128},
    {Sign::NEG, -145, 0xbfffdc00'08fffd78'00c26629'a679ed3b_u128},
    {Sign::NEG, -145, 0xafffe1c0'06eea8e1'23287cb9'd3072728_u128},
    {Sign::NEG, -145, 0x9fffe700'0535541c'd5a37540'fd057315_u128},
    {Sign::NEG, -145, 0x8fffebc0'03cbff32'f82e21c1'fce36810_u128},
    {Sign::NEG, -146, 0xffffe000'05555455'5588887d'dde02702_u128},
    {Sign::NEG, -146, 0xdfffe780'0392aa14'9ac4ed72'adf5b295_u128},
    {Sign::NEG, -146, 0xbfffee00'023fffaf'000c2664'8066b482_u128},
    {Sign::NEG, -146, 0x9ffff380'014d552e'455a3754'b292c077_u128},
    {Sign::NEG, -147, 0xfffff000'01555535'55588888'33333c58_u128},
    {Sign::NEG, -147, 0xbffff700'008ffff5'e000c266'5736679f_u128},
    {Sign::NEG, -148, 0xfffff800'00555551'55558888'85ddde02_u128},
    {Sign::NEG, -149, 0xfffffc00'00155554'd5555888'88733334_u128},
    {Sign::POS, 0, 0_u128},
    {Sign::POS, -148, 0x80000200'000aaaaa'eaaaac44'444eeeef_u128},
    {Sign::POS, -147, 0x80000400'002aaaac'aaaac444'459999ac_u128},
    {Sign::POS, -147, 0xc0000900'0090000a'2000c266'7596679f_u128},
    {Sign::POS, -146, 0x80000800'00aaaaba'aaac4444'6eeef381_u128},
    {Sign::POS, -146, 0xa0000c80'014d557c'655a3755'f81815cc_u128},
    {Sign::POS, -146, 0xc0001200'02400051'000c2668'4c66b482_u128},
    {Sign::POS, -146, 0xe0001880'0392ab40'bac4ed7c'40fb07eb_u128},
    {Sign::POS, -145, 0x80001000'02aaab2a'aac44449'999abe2c_u128},
    {Sign::POS, -145, 0x90001440'03cc00cd'082e21d7'9cbb6812_u128},
    {Sign::POS, -145, 0xa0001900'0535568d'd5a37569'adb01dc3_u128},
    {Sign::POS, -145, 0xb0001e40'06eeac74'33287d01'e8c9d1d9_u128},
    {Sign::POS, -145, 0xc0002400'09000288'00c266a3'2679ed48_u128},
    {Sign::POS, -145, 0xd0002a40'0b7158d1'de776851'22b2764b_u128},
    {Sign::POS, -145, 0xe0003100'0e4aaf5b'2c4ed810'a8063f03_u128},
    {Sign::POS, -145, 0xf0003840'1194062e'0a5143e7'be891c8f_u128},
    {Sign::POS, -144, 0x80002000'0aaaaeaa'ac4444ee'ef3813a1_u128},
    {Sign::POS, -144, 0x88002420'0ccb5a6e'5b7ff7fe'1339025b_u128},
    {Sign::POS, -144, 0x90002880'0f300668'42e21e26'caf39e33_u128},
    {Sign::POS, -144, 0x98002d20'11dcb29e'f271e66f'a5554bc6_u128},
    {Sign::POS, -144, 0xa0003200'14d55f19'5a3757e0'615cc676_u128},
    {Sign::POS, -144, 0xa8003720'181e0bde'ca3b5d82'10ca5cab_u128},
    {Sign::POS, -144, 0xb0003c80'1bbab8f6'f287d25f'3cb032bb_u128},
    {Sign::POS, -144, 0xb8004220'1faf6669'e3278d84'0be28cdb_u128},
    {Sign::POS, -144, 0xc0004800'24001440'0c266dfe'6b482076_u128},
    {Sign::POS, -144, 0xc8004e20'28b0c282'3d9166de'380a6d3d_u128},
    {Sign::POS, -144, 0xd0005480'2dc57139'a7768b35'6ba61e4b_u128},
    {Sign::POS, -144, 0xd8005b20'3342206f'd9e51a18'49db73c1_u128},
    {Sign::POS, -144, 0xe0006200'392ad02e'c4ed8a9d'907eb521_u128},
    {Sign::POS, -144, 0xe8006920'3f838080'b8a197de'a928acd7_u128},
    {Sign::POS, -144, 0xf0007080'46503170'65144cf7'dcc72d3b_u128},
    {Sign::POS, -144, 0xf8007820'4d94e308'da5a1108'890d9f6a_u128},
    {Sign::POS, -143, 0x80004000'2aaacaaa'c4445999'abe2ce2c_u128},
    {Sign::POS, -143, 0x84004410'2ecb2431'1fdbbb4f'3bffc832_u128},
    {Sign::POS, -143, 0x88004840'332d7e1d'97ff8f39'ec91b4ee_u128},
    {Sign::POS, -143, 0x8c004c90'37d3d876'74bcfcf0'b3f0a95d_u128},
    {Sign::POS, -143, 0x90005100'3cc03342'2e21f80c'a6813aff_u128},
    {Sign::POS, -143, 0x94005590'41f48e87'6c3d4629'170ce87f_u128},
    {Sign::POS, -143, 0x98005a40'4772ea4d'071e84e3'b80a8881_u128},
    {Sign::POS, -143, 0x9c005f10'4d3d469a'06d62fdc'bdd6bec3_u128},
    {Sign::POS, -143, 0xa0006400'5355a375'a375a6b7'01dc77c0_u128},
    {Sign::POS, -143, 0xa4006910'59be00e7'450f3318'26ad6b05_u128},
    {Sign::POS, -143, 0xa8006e40'60785ef6'83b60ea8'bd0aa459_u128},
    {Sign::POS, -143, 0xac007390'6786bdab'277e6914'69dd13f5_u128},
    {Sign::POS, -143, 0xb0007900'6eeb1d0d'287d6e0a'0d1e25eb_u128},
    {Sign::POS, -143, 0xb4007e90'76a77d24'aec94b3b'e9b060f5_u128},
    {Sign::POS, -143, 0xb8008440'7ebdddfa'1279365f'ce280cce_u128},
    {Sign::POS, -143, 0xbc008a10'87303f95'dba5732f'3e83e04a_u128},
    {Sign::POS, -143, 0xc0009000'9000a200'c2675967'9ed5b754_u128},
    {Sign::POS, -143, 0xc4009610'99310543'aed95aca'5edb5109_u128},
    {Sign::POS, -143, 0xc8009c40'a2c36967'b917091d'2687160f_u128},
    {Sign::POS, -143, 0xcc00a290'acb9ce76'293d1c2a'0378e75d_u128},
    {Sign::POS, -143, 0xd000a900'b7163478'776977bf'9766f5a7_u128},
    {Sign::POS, -143, 0xd400af90'c1da9b78'4bbb31b1'4776a18b_u128},
    {Sign::POS, -143, 0xd800b640'cd09037f'7e5297d7'6c8564ba_u128},
    {Sign::POS, -143, 0xdc00bd10'd8a36c98'1751360f'8461c447_u128},
    {Sign::POS, -143, 0xe000c400'e4abd6cc'4ed9dc3c'63f44c41_u128},
    {Sign::POS, -143, 0xe400cb10'f1244226'8d10a446'6a5894d5_u128},
    {Sign::POS, -143, 0xe800d240'fe0eaeb1'6a1af81b'b4e6510e_u128},
    {Sign::POS, -143, 0xec00d991'0b6d1c77'ae1f97b0'542a677a_u128},
    {Sign::POS, -143, 0xf000e101'19418b84'51469efe'81d014cc_u128},
    {Sign::POS, -143, 0xf400e891'278dfbe2'7bb98c06'd77a18b4_u128},
    {Sign::POS, -143, 0xf800f041'36546d9d'85a344d0'868bed17_u128},
    {Sign::POS, -143, 0xfc00f811'4596e0c0'f7301d69'90e307cc_u128},
    {Sign::POS, -142, 0x80008000'aaabaaac'4446eef3'8140138f_u128},
    {Sign::POS, -142, 0x82008408'b2cbe5b8'10f5e432'96105497_u128},
    {Sign::POS, -142, 0x84008820'bb2d2189'edbd4f83'ef63f730_u128},
    {Sign::POS, -142, 0x86008c48'c3d05e27'feb654fd'541c638e_u128},
    {Sign::POS, -142, 0x88009080'ccb69b98'7ffadeb8'882f7674_u128},
    {Sign::POS, -142, 0x8a0094c8'd5e0d9e1'c5a59fd3'6bd44397_u128},
};

// 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
constexpr Float128 BIG_COEFFS[4]{
    {Sign::POS, -130, 0xccccccd7'4818e397'7ed78465'd460315b_u128},
    {Sign::NEG, -129, 0x80000000'000478b0'c6388a23'871ce156_u128},
    {Sign::POS, -129, 0xaaaaaaaa'aaaaaaaa'aa807bd8'67763262_u128},
    {Sign::NEG, -128, 0x80000000'00000000'00000000'00000000_u128},
};

[[maybe_unused]] 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);
}
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS

} // namespace

LLVM_LIBC_FUNCTION(double, log1p, (double x)) {
  using FPBits_t = typename fputil::FPBits<double>;

  constexpr int EXP_BIAS = FPBits_t::EXP_BIAS;
  constexpr int FRACTION_LEN = FPBits_t::FRACTION_LEN;
  FPBits_t xbits(x);
  uint64_t x_u = xbits.uintval();

  fputil::DoubleDouble x_dd{0.0, 0.0};

  uint16_t x_exp = xbits.get_biased_exponent();

  if (x_exp >= EXP_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().uintval())) {
        // 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 FPBits_t::inf(Sign::NEG).get_val();
        }
        if (xbits.is_neg() && !xbits.is_nan()) {
          // x < -1.0
          fputil::set_errno_if_required(EDOM);
          fputil::raise_except_if_required(FE_INVALID);
          return FPBits_t::quiet_nan().get_val();
        }
        // x is +Inf or NaN
        if (xbits.is_inf() && xbits.is_pos())
          return x;

        if (xbits.is_signaling_nan())
          fputil::raise_except_if_required(FE_INVALID);
        return FPBits_t::quiet_nan().get_val();
      }
      x_dd.hi = x;
    } else {
      x_dd = fputil::exact_add(x, 1.0);
    }
  } else {
    // |x| < 1
    if (LIBC_UNLIKELY(xbits.get_biased_exponent() <
                      EXP_BIAS - FRACTION_LEN - 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 (x + x == 0.0)
        return x + x; // Handle FTZ/DAZ correctly.

      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);
  uint64_t xhi_frac = xhi_bits.get_mantissa();
  x_u = xhi_bits.uintval();
  // Range reduction:
  // Find k such that |x_hi - k * 2^-7| <= 2^-8.
  int idx = static_cast<int>((xhi_frac + (1ULL << (FRACTION_LEN - 8))) >>
                             (FRACTION_LEN - 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);

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
  // 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];
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS

  // Scale x_dd by 2^(-xh_bits.get_exponent()).
  int64_t s_u = static_cast<int64_t>(x_u & FPBits_t::EXP_MASK) -
                (static_cast<int64_t>(EXP_BIAS) << FRACTION_LEN);
  // Normalize arguments:
  //   1 <= m_dd.hi < 2
  //   |m_dd.lo| < 2^-52.
  // This is exact.
  uint64_t m_hi = FPBits_t::one().uintval() | xhi_frac;

  uint64_t m_lo =
      FPBits_t(x_dd.lo).abs().get_val() > x_dd.hi * 0x1.0p-127
          ? static_cast<uint64_t>(cpp::bit_cast<int64_t>(x_dd.lo) - s_u)
          : 0;

  fputil::DoubleDouble m_dd{FPBits_t(m_lo).get_val(), FPBits_t(m_hi).get_val()};

  // 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_DOUBLE
  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) << (FRACTION_LEN - 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_DOUBLE

  // 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);

#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
  return r1.hi + p;
#else
  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);
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}

} // namespace LIBC_NAMESPACE_DECL