This is an archive of the discontinued LLVM Phabricator instance.

Differential D151049

[libc][math] Implement double precision log1p correctly rounded to all rounding modes.
ClosedPublic

Authored by lntue on May 21 2023, 11:12 AM.

Details

Reviewers
michaelrj
sivachandra
renyichen
santoshn
zimmermann6
cqlauter
Commits
rGb91e78da3780: [libc][math] Implement double precision log1p correctly rounded to all rounding…
Summary

Implement double precision log1p function correctly rounded to all
rounding modes.

Performance

  • For 0.5 <= x <= 2, the fast pass hitting rate is about 99.93%.
  • Benchmarks with ./perf.sh tool from the CORE-MATH project, unit is (CPU clocks / call).
  • Reciprocal throughput from CORE-MATH's perf tool on Ryzen 5900X:
$ ./perf.sh log1p
GNU libc version: 2.35
GNU libc release: stable

-- CORE-MATH reciprocal throughput -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 39.792 + 1.011 clc/call; Median-Min = 0.940 clc/call; Max = 41.373 clc/call;

-- CORE-MATH reciprocal throughput -- without FMA (-march=x86-64-v2)
[####################] 100 %
Ntrial = 20 ; Min = 87.285 + 1.135 clc/call; Median-Min = 1.299 clc/call; Max = 89.715 clc/call;

-- System LIBC reciprocal throughput --
[####################] 100 %
Ntrial = 20 ; Min = 20.666 + 0.123 clc/call; Median-Min = 0.125 clc/call; Max = 20.828 clc/call;

-- LIBC reciprocal throughput -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 20.928 + 0.771 clc/call; Median-Min = 0.725 clc/call; Max = 22.767 clc/call;

-- LIBC reciprocal throughput -- without FMA
[####################] 100 %
Ntrial = 20 ; Min = 31.461 + 0.528 clc/call; Median-Min = 0.602 clc/call; Max = 36.809 clc/call;
  • Latency from CORE-MATH's perf tool on Ryzen 5900X:
$ ./perf.sh log1p --latency
GNU libc version: 2.35
GNU libc release: stable

-- CORE-MATH latency -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 77.875 + 0.062 clc/call; Median-Min = 0.051 clc/call; Max = 78.003 clc/call;

-- CORE-MATH latency -- without FMA (-march=x86-64-v2)
[####################] 100 %
Ntrial = 20 ; Min = 101.958 + 1.202 clc/call; Median-Min = 1.325 clc/call; Max = 104.452 clc/call;

-- System LIBC latency --
[####################] 100 %
Ntrial = 20 ; Min = 60.581 + 1.443 clc/call; Median-Min = 1.611 clc/call; Max = 62.285 clc/call;

-- LIBC latency -- with FMA
[####################] 100 %
Ntrial = 20 ; Min = 48.817 + 1.108 clc/call; Median-Min = 1.300 clc/call; Max = 50.282 clc/call;

-- LIBC latency -- without FMA
[####################] 100 %
Ntrial = 20 ; Min = 61.121 + 0.599 clc/call; Median-Min = 0.761 clc/call; Max = 62.020 clc/call;
  • Accurate pass latency:
$ ./perf.sh log1p --latency --simple_stat
GNU libc version: 2.35
GNU libc release: stable

-- CORE-MATH latency -- with FMA
760.444

-- CORE-MATH latency -- without FMA (-march=x86-64-v2)
827.880

-- LIBC latency -- with FMA
711.837

-- LIBC latency -- without FMA
764.317

Diff Detail

Event Timeline

lntue created this revision.May 21 2023, 11:12 AM
Herald added projects: Restricted Project, Restricted Project. · View Herald TranscriptMay 21 2023, 11:12 AM
Herald added subscribers: ecnelises, tschuett. · View Herald Transcript
lntue requested review of this revision.May 21 2023, 11:12 AM
Harbormaster completed remote builds in B233441: Diff 524117.May 21 2023, 11:12 AM
lntue edited the summary of this revision. (Show Details)May 21 2023, 11:44 AM
Herald added a subscriber: pengfei. · View Herald TranscriptMay 21 2023, 11:44 AM
lntue added a parent revision: D150374: [libc][math] Implement double precision log2 function correctly rounded to all rounding modes..May 21 2023, 11:44 AM
lntue updated this revision to Diff 524122.May 21 2023, 11:44 AM

Sync.

Harbormaster completed remote builds in B233444: Diff 524122.May 21 2023, 11:50 AM
lntue updated this revision to Diff 524150.May 21 2023, 8:00 PM

Sync.

Harbormaster completed remote builds in B233465: Diff 524150.May 21 2023, 8:07 PM
zimmermann6 requested changes to this revision.May 23 2023, 6:35 AM

I get the same error as for log2:

CMake Error at /localdisk/zimmerma/llvm-project/libc/cmake/modules/LLVMLibCLibraryRules.cmake:5 (get_target_property):
  get_target_property() called with non-existent target
  "libc.src.math.generic.log_range_reduction".
Call Stack (most recent call first):
  /localdisk/zimmerma/llvm-project/libc/cmake/modules/LLVMLibCLibraryRules.cmake:35 (collect_object_file_deps)
  /localdisk/zimmerma/llvm-project/libc/cmake/modules/LLVMLibCLibraryRules.cmake:35 (collect_object_file_deps)
  /localdisk/zimmerma/llvm-project/libc/cmake/modules/LLVMLibCLibraryRules.cmake:82 (collect_object_file_deps)
  /localdisk/zimmerma/llvm-project/libc/lib/CMakeLists.txt:26 (add_entrypoint_library)
This revision now requires changes to proceed.May 23 2023, 6:35 AM
lntue updated this revision to Diff 524717.May 23 2023, 7:50 AM

Sync to HEAD.

zimmermann6 accepted this revision.May 23 2023, 8:01 AM

this function is faster than core-math, even for the reciprocal throughput, great work!

GNU libc version: 2.36
GNU libc release: stable
[####################] 100 %
Ntrial = 20 ; Min = 49.596 + 0.352 clc/call; Median-Min = 0.315 clc/call; Max = 50.216 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 26.743 + 0.357 clc/call; Median-Min = 0.333 clc/call; Max = 29.183 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 38.887 + 0.296 clc/call; Median-Min = 0.249 clc/call; Max = 41.140 clc/call;
GNU libc version: 2.36
GNU libc release: stable
[####################] 100 %
Ntrial = 20 ; Min = 94.356 + 0.365 clc/call; Median-Min = 0.288 clc/call; Max = 95.270 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 69.370 + 0.356 clc/call; Median-Min = 0.262 clc/call; Max = 70.010 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 67.833 + 0.346 clc/call; Median-Min = 0.292 clc/call; Max = 68.525 clc/call;
This revision is now accepted and ready to land.May 23 2023, 8:01 AM
This revision was landed with ongoing or failed builds.May 23 2023, 8:04 AM
Closed by commit rGb91e78da3780: [libc][math] Implement double precision log1p correctly rounded to all rounding… (authored by lntue). · Explain Why
This revision was automatically updated to reflect the committed changes.
lntue added a commit: rGb91e78da3780: [libc][math] Implement double precision log1p correctly rounded to all rounding….
Harbormaster completed remote builds in B233872: Diff 524717.May 23 2023, 8:18 AM

Revision Contents

PathSize
libc/
config/
darwin/
arm/
1 line
linux/
aarch64/
1 line
x86_64/
1 line
windows/
1 line
spec/
1 line
src/
math/
1 line
generic/
20 lines
1038 lines
18 lines
test/
src/
math/
14 lines
166 lines

Diff 524726

libc/config/darwin/arm/entrypoints.txt

Show First 20 Lines • Show All 170 Lines • ▼ Show 20 Linesset(TARGET_LIBM_ENTRYPOINTS
libc.src.math.lround libc.src.math.lround
libc.src.math.lroundf libc.src.math.lroundf
libc.src.math.lroundl libc.src.math.lroundl
libc.src.math.ldexp libc.src.math.ldexp
libc.src.math.ldexpf libc.src.math.ldexpf
libc.src.math.ldexpl libc.src.math.ldexpl
libc.src.math.log10 libc.src.math.log10
libc.src.math.log10f libc.src.math.log10f
libc.src.math.log1p
libc.src.math.log1pf libc.src.math.log1pf
libc.src.math.log2 libc.src.math.log2
libc.src.math.log2f libc.src.math.log2f
libc.src.math.log libc.src.math.log
libc.src.math.logf libc.src.math.logf
libc.src.math.logb libc.src.math.logb
libc.src.math.logbf libc.src.math.logbf
libc.src.math.logbl libc.src.math.logbl
▲ Show 20 Lines • Show All 41 Lines • Show Last 20 Lines

libc/config/linux/aarch64/entrypoints.txt

Show First 20 Lines • Show All 281 Lines • ▼ Show 20 Linesset(TARGET_LIBM_ENTRYPOINTS
libc.src.math.lround libc.src.math.lround
libc.src.math.lroundf libc.src.math.lroundf
libc.src.math.lroundl libc.src.math.lroundl
libc.src.math.ldexp libc.src.math.ldexp
libc.src.math.ldexpf libc.src.math.ldexpf
libc.src.math.ldexpl libc.src.math.ldexpl
libc.src.math.log10 libc.src.math.log10
libc.src.math.log10f libc.src.math.log10f
libc.src.math.log1p
libc.src.math.log1pf libc.src.math.log1pf
libc.src.math.log2 libc.src.math.log2
libc.src.math.log2f libc.src.math.log2f
libc.src.math.log libc.src.math.log
libc.src.math.logf libc.src.math.logf
libc.src.math.logb libc.src.math.logb
libc.src.math.logbf libc.src.math.logbf
libc.src.math.logbl libc.src.math.logbl
▲ Show 20 Lines • Show All 190 Lines • Show Last 20 Lines

libc/config/linux/x86_64/entrypoints.txt

Show First 20 Lines • Show All 286 Lines • ▼ Show 20 Linesset(TARGET_LIBM_ENTRYPOINTS
libc.src.math.llrint libc.src.math.llrint
libc.src.math.llrintf libc.src.math.llrintf
libc.src.math.llrintl libc.src.math.llrintl
libc.src.math.llround libc.src.math.llround
libc.src.math.llroundf libc.src.math.llroundf
libc.src.math.llroundl libc.src.math.llroundl
libc.src.math.log10 libc.src.math.log10
libc.src.math.log10f libc.src.math.log10f
libc.src.math.log1p
libc.src.math.log1pf libc.src.math.log1pf
libc.src.math.log2 libc.src.math.log2
libc.src.math.log2f libc.src.math.log2f
libc.src.math.log libc.src.math.log
libc.src.math.logf libc.src.math.logf
libc.src.math.logb libc.src.math.logb
libc.src.math.logbf libc.src.math.logbf
libc.src.math.logbl libc.src.math.logbl
▲ Show 20 Lines • Show All 231 Lines • Show Last 20 Lines

libc/config/windows/entrypoints.txt

Show First 20 Lines • Show All 163 Lines • ▼ Show 20 Linesset(TARGET_LIBM_ENTRYPOINTS
libc.src.math.llrint libc.src.math.llrint
libc.src.math.llrintf libc.src.math.llrintf
libc.src.math.llrintl libc.src.math.llrintl
libc.src.math.llround libc.src.math.llround
libc.src.math.llroundf libc.src.math.llroundf
libc.src.math.llroundl libc.src.math.llroundl
libc.src.math.log10 libc.src.math.log10
libc.src.math.log10f libc.src.math.log10f
libc.src.math.log1p
libc.src.math.log1pf libc.src.math.log1pf
libc.src.math.log2 libc.src.math.log2
libc.src.math.log2f libc.src.math.log2f
libc.src.math.log libc.src.math.log
libc.src.math.logf libc.src.math.logf
libc.src.math.logb libc.src.math.logb
libc.src.math.logbf libc.src.math.logbf
libc.src.math.logbl libc.src.math.logbl
▲ Show 20 Lines • Show All 50 Lines • Show Last 20 Lines

libc/spec/stdc.td

Show First 20 Lines • Show All 402 Lines • ▼ Show 20 Lines HeaderSpec Math = HeaderSpec<
FunctionSpec<"ldexp", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<IntType>]>, FunctionSpec<"ldexp", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<IntType>]>,
FunctionSpec<"ldexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntType>]>, FunctionSpec<"ldexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntType>]>,
FunctionSpec<"ldexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<IntType>]>, FunctionSpec<"ldexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<IntType>]>,
FunctionSpec<"log10", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>, FunctionSpec<"log10", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"log10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>, FunctionSpec<"log10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"log1p", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"log1pf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>, FunctionSpec<"log1pf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"log2", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>, FunctionSpec<"log2", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"log2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>, FunctionSpec<"log2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"log", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>, FunctionSpec<"log", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"logf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>, FunctionSpec<"logf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
▲ Show 20 Lines • Show All 698 Lines • Show Last 20 Lines

libc/src/math/CMakeLists.txt

Show First 20 Lines • Show All 108 Lines • ▼ Show 20 Lines
add_math_entrypoint_object(ldexp) add_math_entrypoint_object(ldexp)
add_math_entrypoint_object(ldexpf) add_math_entrypoint_object(ldexpf)
add_math_entrypoint_object(ldexpl) add_math_entrypoint_object(ldexpl)
add_math_entrypoint_object(log10) add_math_entrypoint_object(log10)
add_math_entrypoint_object(log10f) add_math_entrypoint_object(log10f)
add_math_entrypoint_object(log1p)
add_math_entrypoint_object(log1pf) add_math_entrypoint_object(log1pf)
add_math_entrypoint_object(log2) add_math_entrypoint_object(log2)
add_math_entrypoint_object(log2f) add_math_entrypoint_object(log2f)
add_math_entrypoint_object(log) add_math_entrypoint_object(log)
add_math_entrypoint_object(logf) add_math_entrypoint_object(logf)
▲ Show 20 Lines • Show All 69 Lines • Show Last 20 Lines

libc/src/math/generic/CMakeLists.txt

Show First 20 Lines • Show All 809 Lines • ▼ Show 20 Lines DEPENDS
libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.fma libc.src.__support.FPUtil.fma
libc.src.__support.FPUtil.polyeval libc.src.__support.FPUtil.polyeval
COMPILE_OPTIONS COMPILE_OPTIONS
-O3 -O3
) )
add_entrypoint_object( 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 log1pf
SRCS SRCS
log1pf.cpp log1pf.cpp
HDRS HDRS
../log1pf.h ../log1pf.h
DEPENDS DEPENDS
.common_constants .common_constants
libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.except_value_utils
▲ Show 20 Lines • Show All 705 Lines • Show Last 20 Lines

libc/src/math/generic/log1p.cpp

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

libc/src/math/log1p.h

  • This file was added.
//===-- 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

libc/test/src/math/CMakeLists.txt

Show First 20 Lines • Show All 1,346 Lines • ▼ Show 20 Linesadd_fp_unittest(
DEPENDS DEPENDS
libc.src.errno.errno libc.src.errno.errno
libc.include.math libc.include.math
libc.src.math.log10f libc.src.math.log10f
libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.fp_bits
) )
add_fp_unittest( 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 log1pf_test
NEED_MPFR NEED_MPFR
SUITE SUITE
libc_math_unittests libc_math_unittests
SRCS SRCS
log1pf_test.cpp log1pf_test.cpp
DEPENDS DEPENDS
libc.src.errno.errno libc.src.errno.errno
▲ Show 20 Lines • Show All 242 Lines • Show Last 20 Lines

libc/test/src/math/log1p_test.cpp

  • This file was added.
//===-- Unittests for log1p -----------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#include "src/__support/FPUtil/FPBits.h"
#include "src/errno/libc_errno.h"
#include "src/math/log1p.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
#include <math.h>
#include <errno.h>
#include <stdint.h>
namespace mpfr = __llvm_libc::testing::mpfr;
using __llvm_libc::testing::tlog;
DECLARE_SPECIAL_CONSTANTS(double)
TEST(LlvmLibcLog1pTest, SpecialNumbers) {
EXPECT_FP_EQ(aNaN, __llvm_libc::log1p(aNaN));
EXPECT_FP_EQ(inf, __llvm_libc::log1p(inf));
EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log1p(neg_inf), FE_INVALID);
EXPECT_FP_IS_NAN_WITH_EXCEPTION(__llvm_libc::log1p(-2.0), FE_INVALID);
EXPECT_FP_EQ(zero, __llvm_libc::log1p(0.0));
EXPECT_FP_EQ(neg_zero, __llvm_libc::log1p(-0.0));
EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, __llvm_libc::log1p(-1.0), FE_DIVBYZERO);
}
TEST(LlvmLibcLog1pTest, TrickyInputs) {
constexpr int N = 41;
constexpr uint64_t INPUTS[N] = {
0x3ff0000000000000, // x = 1.0
0x4024000000000000, // x = 10.0
0x4059000000000000, // x = 10^2
0x408f400000000000, // x = 10^3
0x40c3880000000000, // x = 10^4
0x40f86a0000000000, // x = 10^5
0x412e848000000000, // x = 10^6
0x416312d000000000, // x = 10^7
0x4197d78400000000, // x = 10^8
0x41cdcd6500000000, // x = 10^9
0x4202a05f20000000, // x = 10^10
0x42374876e8000000, // x = 10^11
0x426d1a94a2000000, // x = 10^12
0x42a2309ce5400000, // x = 10^13
0x42d6bcc41e900000, // x = 10^14
0x430c6bf526340000, // x = 10^15
0x4341c37937e08000, // x = 10^16
0x4376345785d8a000, // x = 10^17
0x43abc16d674ec800, // x = 10^18
0x43e158e460913d00, // x = 10^19
0x4415af1d78b58c40, // x = 10^20
0x444b1ae4d6e2ef50, // x = 10^21
0x4480f0cf064dd592, // x = 10^22
0x3fefffffffef06ad, 0x3fefde0f22c7d0eb, 0x225e7812faadb32f,
0x3fee1076964c2903, 0x3fdfe93fff7fceb0, 0x3ff012631ad8df10,
0x3fefbfdaa448ed98, 0x3fd00a8cefe9a5f8, 0x3fd0b4d870eb22f8,
0x3c90c40cef04efb5, 0x449d2ccad399848e, 0x4aa12ccdffd9d2ec,
0x5656f070b92d36ce, 0x6db06dcb74f76bcc, 0x7f1954e72ffd4596,
0x5671e2f1628093e4, 0x73dac56e2bf1a951, 0x8001bc6879ea14c5,
};
for (int i = 0; i < N; ++i) {
double x = double(FPBits(INPUTS[i]));
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log1p, x,
__llvm_libc::log1p(x), 0.5);
}
}
TEST(LlvmLibcLog1pTest, AllExponents) {
double x = 0x1.0p-1074;
for (int i = -1074; i < 1024; ++i, x *= 2.0) {
ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Log1p, x,
__llvm_libc::log1p(x), 0.5);
}
}
TEST(LlvmLibcLog1pTest, InDoubleRange) {
constexpr uint64_t COUNT = 1234561;
auto test = [&](uint64_t start, uint64_t stop,
mpfr::RoundingMode rounding_mode) {
mpfr::ForceRoundingMode __r(rounding_mode);
uint64_t fails = 0;
uint64_t count = 0;
uint64_t cc = 0;
double mx, mr = 0.0;
double tol = 0.5;
uint64_t step = (stop - start) / COUNT;
for (uint64_t i = 0, v = start; i <= COUNT; ++i, v += step) {
double x = FPBits(v).get_val();
if (isnan(x) || isinf(x) || x < 0.0)
continue;
libc_errno = 0;
double result = __llvm_libc::log1p(x);
++cc;
if (isnan(result) || isinf(result))
continue;
++count;
// ASSERT_MPFR_MATCH(mpfr::Operation::Log1p, x, result, 0.5);
if (!EXPECT_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Log1p, x,
result, 0.5, rounding_mode)) {
++fails;
while (!EXPECT_MPFR_MATCH_ROUNDING_SILENTLY(
mpfr::Operation::Log1p, x, result, tol, rounding_mode)) {
mx = x;
mr = result;
tol *= 2.0;
}
}
}
tlog << " Log1p failed: " << fails << "/" << count << "/" << cc
<< " tests.\n";
tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
if (fails) {
EXPECT_MPFR_MATCH(mpfr::Operation::Log1p, mx, mr, 0.5, rounding_mode);
}
};
auto test_all_rounding = [&](uint64_t start, uint64_t stop,
const char *start_str, const char *stop_str) {
tlog << "\n=== Test in range [" << start_str << ", " << stop_str
<< "] ===\n";
tlog << "\n Test Rounding To Nearest...\n";
test(start, stop, mpfr::RoundingMode::Nearest);
tlog << "\n Test Rounding Downward...\n";
test(start, stop, mpfr::RoundingMode::Downward);
tlog << "\n Test Rounding Upward...\n";
test(start, stop, mpfr::RoundingMode::Upward);
tlog << "\n Test Rounding Toward Zero...\n";
test(start, stop, mpfr::RoundingMode::TowardZero);
};
test_all_rounding(0x0000'0000'0000'0001ULL, 0x0010'0000'0000'0000ULL,
"2^-1074", "2^-1022");
test_all_rounding(0x39B0'0000'0000'0000ULL, 0x3A50'0000'0000'0000ULL,
"2^-100", "2^-90");
test_all_rounding(0x3CD0'0000'0000'0000ULL, 0x3D20'0000'0000'0000ULL, "2^-50",
"2^-45");
test_all_rounding(0x3E10'0000'0000'0000ULL, 0x3E40'0000'0000'0000ULL, "2^-30",
"2^-27");
test_all_rounding(0x3FD0'0000'0000'0000ULL, 0x4010'0000'0000'0000ULL, "0.25",
"4.0");
test_all_rounding(0x4630'0000'0000'0000ULL, 0x4670'0000'0000'0000ULL, "2^100",
"2^104");
test_all_rounding(0x7FD0'0000'0000'0000ULL, 0x7FF0'0000'0000'0000ULL,
"2^1022", "2^1024");
}