llvm/libc/src/math/generic/range_reduction_double_common.h

//===-- Range reduction for double precision sin/cos/tan -*- 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_GENERIC_RANGE_REDUCTION_DOUBLE_COMMON_H
#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_COMMON_H

#include "src/__support/FPUtil/FPBits.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/common.h"
#include "src/__support/integer_literals.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"

namespace LIBC_NAMESPACE_DECL {

#ifdef LIBC_TARGET_CPU_HAS_FMA
static constexpr unsigned SPLIT = DEFAULT_DOUBLE_SPLIT;
#else
// When there is no-FMA instructions, in order to have exact product of 2 double
// precision with directional roundings, we need to lower the precision of the
// constants by at least 1 bit, and use a different splitting constant.
static constexpr unsigned SPLIT = …;
#endif // LIBC_TARGET_CPU_HAS_FMA

DoubleDouble;
Float128;

#define FAST_PASS_EXPONENT …

// For 2^-7 < |x| < 2^16, return k and u such that:
//   k = round(x * 128/pi)
//   x mod pi/128 = x - k * pi/128 ~ u.hi + u.lo
// Error bound:
//   |(x - k * pi/128) - (u_hi + u_lo)| <= max(ulp(ulp(u_hi)), 2^-119)
//                                      <= 2^-111.
LIBC_INLINE unsigned range_reduction_small(double x, DoubleDouble &u) { … }

// Digits of 2^(16*i) / pi, generated by Sollya with:
// > procedure ulp(x, n) { return 2^(floor(log2(abs(x))) - n); };
// > for i from 0 to 63 do {
//     if i < 3 then { pi_inv = 0.25 + 2^(16*(i - 3)) / pi; }
//     else { pi_inv = 2^(16*(i-3)) / pi; };
//     pn = nearestint(pi_inv);
//     pi_frac = pi_inv - pn;
//     a = round(pi_frac, 51, RN);
//     b = round(pi_frac - a, 51, RN);
//     c = round(pi_frac - a - b, 51, RN);
//     d = round(pi_frac - a - b - c, D, RN);
//     print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},");
// };
//
// Notice that for [0..2] the leading bit of 2^(16*(i - 3)) / pi is very small,
// so we add 0.25 so that the conditions for the algorithms are still satisfied,
// and one of those conditions guarantees that ulp(0.25 * x_reduced) >= 2, and
// will safely be discarded.

static constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = …;

// For large range |x| >= 2^16, we perform the range reduction computations as:
//   u = x - k * pi/128 = (pi/128) * (x * (128/pi) - k).
// We use the exponent of x to find 4 double-chunks of 128/pi:
// c_hi, c_mid, c_lo, c_lo_2 such that:
//   1) ulp(round(x * c_hi, D, RN)) >= 2^8 = 256,
//   2) If x * c_hi = ph_hi + ph_lo and x * c_mid = pm_hi + pm_lo, then
//        min(ulp(ph_lo), ulp(pm_hi)) >= 2^-53.
// This will allow us to drop the high part ph_hi and the addition:
//   (ph_lo + pm_hi) mod 1
// can be exactly representable in a double precision.
// This will allow us to do split the computations as:
//   (x * 256/pi) ~ x * (c_hi + c_mid + c_lo + c_lo_2)    (mod 256)
//                ~ (ph_lo + pm_hi) + (pm_lo + x * c_lo) + x * c_lo_2.
// Then,
//   round(x * 128/pi) = round(ph_lo + pm_hi)    (mod 256)
// And the high part of fractional part of (x * 128/pi) can simply be:
//   {x * 128/pi}_hi = {ph_lo + pm_hi}.
// To prevent overflow when x is very large, we simply scale up
// (c_hi, c_mid, c_lo, c_lo_2) by a fixed power of 2 (based on the index) and
// scale down x by the same amount.

struct LargeRangeReduction { … };

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
static Float128 range_reduction_small_f128(double x) { … }

static constexpr Float128 SIN_K_PI_OVER_128_F128[65] = …;
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_COMMON_H