log1p.cpp | Explore in Territory

//===-- 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 "common_constants.h"

namespace LIBC_NAMESPACE_DECL {

// 128-bit precision dyadic floating point numbers.
Float128;

operator""_u128;

namespace {

// Extra errors from P is from using x^2 to reduce evaluation latency and
// directional rounding.
constexpr double P_ERR = …;

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

// R1[i] = 2^-8 * nearestint( 2^8 / (1 + i * 2^-7) )
constexpr double R1[129] = …;

// 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] = …;

// 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] = …;

// 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] = …;

// -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] = …;

// 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] = …;

// -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] = …;

// 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] = …;

// -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] = …;

// 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]{ … };

LIBC_INLINE double log1p_accurate(int e_x, int index,
                                  fputil::DoubleDouble m_x) { … }

} // namespace

LLVM_LIBC_FUNCTION(double, log1p, (double x)) { … }

} // namespace LIBC_NAMESPACE_DECL
llvm/libc/src/math/generic/log1p.cpp
