expm1.cpp | Explore in Territory

//===-- Double-precision e^x - 1 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/expm1.h"
#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
#include "explogxf.h"         // ziv_test_denorm.
#include "src/__support/CPP/bit.h"
#include "src/__support/CPP/optional.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/except_value_utils.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/FPUtil/triple_double.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

#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0)
#define LIBC_MATH_EXPM1_SKIP_ACCURATE_PASS
#endif

namespace LIBC_NAMESPACE_DECL {

DoubleDouble;
TripleDouble;
Float128;

operator""_u128;

// log2(e)
constexpr double LOG2_E = …;

// Error bounds:
// Errors when using double precision.
// 0x1.8p-63;
constexpr uint64_t ERR_D = …;
// Errors when using double-double precision.
// 0x1.0p-99
[[maybe_unused]] constexpr uint64_t ERR_DD = …;

// -2^-12 * log(2)
// > a = -2^-12 * log(2);
// > b = round(a, 30, RN);
// > c = round(a - b, 30, RN);
// > d = round(a - b - c, D, RN);
// Errors < 1.5 * 2^-133
constexpr double MLOG_2_EXP2_M12_HI = …;
constexpr double MLOG_2_EXP2_M12_MID = …;
constexpr double MLOG_2_EXP2_M12_MID_30 = …;
constexpr double MLOG_2_EXP2_M12_LO = …;

namespace {

// Polynomial approximations with double precision:
// Return expm1(dx) / x ~ 1 + dx / 2 + dx^2 / 6 + dx^3 / 24.
// For |dx| < 2^-13 + 2^-30:
//   | output - expm1(dx) / dx | < 2^-51.
LIBC_INLINE double poly_approx_d(double dx) { … }

// Polynomial approximation with double-double precision:
// Return expm1(dx) / dx ~ 1 + dx / 2 + dx^2 / 6 + ... + dx^6 / 5040
// For |dx| < 2^-13 + 2^-30:
//   | output - expm1(dx) | < 2^-101
DoubleDouble poly_approx_dd(const DoubleDouble &dx) { … }

// Polynomial approximation with 128-bit precision:
// Return (exp(dx) - 1)/dx ~ 1 + dx / 2 + dx^2 / 6 + ... + dx^6 / 5040
// For |dx| < 2^-13 + 2^-30:
//   | output - exp(dx) | < 2^-126.
[[maybe_unused]] Float128 poly_approx_f128(const Float128 &dx) { … }

#ifdef DEBUGDEBUG
std::ostream &operator<<(std::ostream &OS, const Float128 &r) {
  OS << (r.sign == Sign::NEG ? "-(" : "(") << r.mantissa.val[0] << " + "
     << r.mantissa.val[1] << " * 2^64) * 2^" << r.exponent << "\n";
  return OS;
}

std::ostream &operator<<(std::ostream &OS, const DoubleDouble &r) {
  OS << std::hexfloat << "(" << r.hi << " + " << r.lo << ")"
     << std::defaultfloat << "\n";
  return OS;
}
#endif

// Compute exp(x) - 1 using 128-bit precision.
// TODO(lntue): investigate triple-double precision implementation for this
// step.
[[maybe_unused]] Float128 expm1_f128(double x, double kd, int idx1, int idx2) { … }

// Compute exp(x) - 1 with double-double precision.
DoubleDouble exp_double_double(double x, double kd, const DoubleDouble &exp_mid,
                               const DoubleDouble &hi_part) { … }

// Check for exceptional cases when
// |x| <= 2^-53 or x < log(2^-54) or x >= 0x1.6232bdd7abcd3p+9
double set_exceptional(double x) { … }

} // namespace

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

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