llvm/libc/src/math/generic/range_reduction_double_nofma.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_NOFMA_H
#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_NOFMA_H

#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"

namespace LIBC_NAMESPACE_DECL {

namespace nofma {

DoubleDouble;

LIBC_INLINE constexpr int FAST_PASS_EXPONENT = …;

// Digits of 2^(16*i) / pi, generated by Sollya with:
// For [2..62]:
// >  for i from 3 to 63 do {
//     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, D, RN);
//     d = round(pi_frac - a - b - c, D, RN);
//     print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},");
//   };
// For [0..1]:
// 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.
//  for i from 0 to 2 do {
//     pi_frac = 0.25 + 2^(16*(i - 3)) / pi;
//     a = round(pi_frac, 51, RN);
//     b = round(pi_frac - a, 51, RN);
//     c = round(pi_frac - a - b, D, RN);
//     d = round(pi_frac - a - b - c, D, RN);
//     print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},");
//   };
// For The fast pass using double-double, we only need 3 parts (a, b, c), but
// for the accurate pass using Float128, instead of using another table of
// Float128s, we simply add the fourth path (a, b, c, d), which simplify the
// implementation a bit and saving some memory.
LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = …;

// Lookup table for sin(k * pi / 128) with k = 0, ..., 255.
// Table is generated with Sollya as follow:
// > display = hexadecimal;
// > for k from 0 to 255 do {
//     a = round(sin(k * pi/128), 51, RN);
//     b = round(sin(k * pi/128) - a, D, RN);
//     print("{", b, ",", a, "},");
//   };
LIBC_INLINE constexpr DoubleDouble SIN_K_PI_OVER_128[256] = …;

LIBC_INLINE unsigned range_reduction_small(double x, DoubleDouble &u) { … }

} // namespace nofma

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_DOUBLE_NOFMA_H