exp.cpp | Explore in Territory

//===-- Double-precision e^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/exp.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/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

#include <errno.h>

namespace LIBC_NAMESPACE_DECL {

DoubleDouble;
TripleDouble;
Float128;

operator""_u128;

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

// Error bounds:
// Errors when using double precision.
constexpr double ERR_D = …;
// Errors when using double-double precision.
constexpr double 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 exp(dx) ~ 1 + dx + dx^2 / 2 + ... + dx^6 / 720
// For |dx| < 2^-13 + 2^-30:
//   | output - exp(dx) | < 2^-101
DoubleDouble poly_approx_dd(const DoubleDouble &dx) { … }

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

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

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

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

} // namespace

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

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