explogxf.h | Explore in Territory

//===-- Single-precision general exp/log functions ------------------------===//
//
// 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_GENERIC_EXPLOGXF_H
#define LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H

#include "common_constants.h"
#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/nearest_integer.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/properties/cpu_features.h"

#include <errno.h>

namespace LIBC_NAMESPACE_DECL {

struct ExpBase { … };

struct Exp10Base : public ExpBase { … };

constexpr int LOG_P1_BITS = …;
constexpr int LOG_P1_SIZE = …;

// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
extern const double LOG_P1_LOG2[LOG_P1_SIZE];

// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
extern const double LOG_P1_1_OVER[LOG_P1_SIZE];

// Taylor series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
// K_LOG2_ODD starts from x^3
extern const double K_LOG2_ODD[4];
extern const double K_LOG2_EVEN[4];

// Output of range reduction for exp_b: (2^(mid + hi), lo)
// where:
//   b^x = 2^(mid + hi) * b^lo
struct exp_b_reduc_t { … };

// The function correctly calculates b^x value with at least float precision
// in a limited range.
// Range reduction:
//   b^x = 2^(hi + mid) * b^lo
// where:
//   x = (hi + mid) * log_b(2) + lo
//   hi is an integer,
//   0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer
//   -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1)
// Base class needs to provide the following constants:
//   - MID_BITS    : number of bits after decimal points used for mid
//   - MID_MASK    : 2^MID_BITS - 1, mask to extract mid bits
//   - LOG2_B      : log2(b) * 2^MID_BITS for scaling
//   - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS)
//   - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS)
//   - EXP_2_MID   : look up table for bit fields of 2^mid
// Return:
//   { 2^(hi + mid), lo }
template <class Base> LIBC_INLINE exp_b_reduc_t exp_b_range_reduc(float x) { … }

// The function correctly calculates sinh(x) and cosh(x) by calculating exp(x)
// and exp(-x) simultaneously.
// To compute e^x, we perform the following range
// reduction: find hi, mid, lo such that:
//   x = (hi + mid) * log(2) + lo, in which
//     hi is an integer,
//     0 <= mid * 2^5 < 32 is an integer
//     -2^(-6) <= lo * log2(e) <= 2^-6.
// In particular,
//   hi + mid = round(x * log2(e) * 2^5) * 2^(-5).
// Then,
//   e^x = 2^(hi + mid) * e^lo = 2^hi * 2^mid * e^lo.
// 2^mid is stored in the lookup table of 32 elements.
// e^lo is computed using a degree-5 minimax polynomial
// generated by Sollya:
//   e^lo ~ P(lo) = 1 + lo + c2 * lo^2 + ... + c5 * lo^5
//        = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)
//        = P_even + lo * P_odd
// We perform 2^hi * 2^mid by simply add hi to the exponent field
// of 2^mid.
// To compute e^(-x), notice that:
//   e^(-x) = 2^(-(hi + mid)) * e^(-lo)
//          ~ 2^(-(hi + mid)) * P(-lo)
//          = 2^(-(hi + mid)) * (P_even - lo * P_odd)
// So:
//   sinh(x) = (e^x - e^(-x)) / 2
//           ~ 0.5 * (2^(hi + mid) * (P_even + lo * P_odd) -
//                    2^(-(hi + mid)) * (P_even - lo * P_odd))
//           = 0.5 * (P_even * (2^(hi + mid) - 2^(-(hi + mid))) +
//                    lo * P_odd * (2^(hi + mid) + 2^(-(hi + mid))))
// And similarly:
//   cosh(x) = (e^x + e^(-x)) / 2
//           ~ 0.5 * (P_even * (2^(hi + mid) + 2^(-(hi + mid))) +
//                    lo * P_odd * (2^(hi + mid) - 2^(-(hi + mid))))
// The main point of these formulas is that the expensive part of calculating
// the polynomials approximating lower parts of e^(x) and e^(-x) are shared
// and only done once.
template <bool is_sinh> LIBC_INLINE double exp_pm_eval(float x) { … }

// x should be positive, normal finite value
LIBC_INLINE static double log2_eval(double x) { … }

// x should be positive, normal finite value
LIBC_INLINE static double log_eval(double x) { … }

// Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We
// assume further that 1 <= mid < 2, mid + lo < 2, and |lo| << mid.
// Notice that, if 0 < x < 2^-1022,
//   double(2^-1022 + x) - 2^-1022 = double(x).
// So if we scale x up by 2^1022, we can use
//   double(1.0 + 2^1022 * x) - 1.0 to test how x is rounded in denormal range.
LIBC_INLINE cpp::optional<double> ziv_test_denorm(int hi, double mid, double lo,
                                                  double err) { … }

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H
llvm/libc/src/math/generic/explogxf.h
