libc/src/math/generic/range_reduction_fma.h

//===-- Utilities for trigonometric functions with FMA ----------*- 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_FMA_H
#define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H

#include "src/__support/FPUtil/FMA.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/common.h"

namespace LIBC_NAMESPACE {

namespace fma {

static constexpr uint32_t FAST_PASS_BOUND = 0x5600'0000U; // 2^45

// Digits of 32/pi, generated by Sollya with:
// > a0 = D(32/pi);
// > a1 = D(32/pi - a0);
// > a2 = D(32/pi - a0 - a1);
// > a3 = D(32/pi - a0 - a1 - a2);
static constexpr double THIRTYTWO_OVER_PI[5] = {
    0x1.45f306dc9c883p+3, -0x1.6b01ec5417056p-51, -0x1.6447e493ad4cep-105,
    0x1.e21c820ff28b2p-159, -0x1.508510ea79237p-214};

// Return k and y, where
//   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
LIBC_INLINE int64_t small_range_reduction(double x, double &y) {
  double kd = fputil::nearest_integer(x * THIRTYTWO_OVER_PI[0]);
  y = fputil::fma(x, THIRTYTWO_OVER_PI[0], -kd);
  y = fputil::fma(x, THIRTYTWO_OVER_PI[1], y);
  return static_cast<int64_t>(kd);
}

// Return k and y, where
//   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
// This is used for sinf, cosf, sincosf.
LIBC_INLINE int64_t large_range_reduction(double x, int x_exp, double &y) {
  // 2^45 <= |x| < 2^99
  if (x_exp < 99) {
    // - When x < 2^99, the full exact product of x * THIRTYTWO_OVER_PI[0]
    // contains at least one integral bit <= 2^5.
    // - When 2^45 <= |x| < 2^55, the lowest 6 unit bits are contained
    // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[0]).
    // - When |x| >= 2^55, the LSB of double(x * THIRTYTWO_OVER_PI[0]) is at
    // least 2^6.
    fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[0]);
    prod_hi.bits &= (x_exp < 55) ? (~0xfffULL) : (~0ULL); // |x| < 2^55
    double k_hi = fputil::nearest_integer(prod_hi.get_val());
    double truncated_prod = fputil::fma(x, THIRTYTWO_OVER_PI[0], -k_hi);
    double prod_lo = fputil::fma(x, THIRTYTWO_OVER_PI[1], truncated_prod);
    double k_lo = fputil::nearest_integer(prod_lo);
    y = fputil::fma(x, THIRTYTWO_OVER_PI[1], truncated_prod - k_lo);
    y = fputil::fma(x, THIRTYTWO_OVER_PI[2], y);
    y = fputil::fma(x, THIRTYTWO_OVER_PI[3], y);

    return static_cast<int64_t>(k_lo);
  }

  // - When x >= 2^110, the full exact product of x * THIRTYTWO_OVER_PI[0] does
  // not contain any of the lowest 6 unit bits, so we can ignore it completely.
  // - When 2^99 <= |x| < 2^110, the lowest 6 unit bits are contained
  // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[1]).
  // - When |x| >= 2^110, the LSB of double(x * THIRTYTWO_OVER_PI[1]) is at
  // least 64.
  fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[1]);
  prod_hi.bits &= (x_exp < 110) ? (~0xfffULL) : (~0ULL); // |x| < 2^110
  double k_hi = fputil::nearest_integer(prod_hi.get_val());
  double truncated_prod = fputil::fma(x, THIRTYTWO_OVER_PI[1], -k_hi);
  double prod_lo = fputil::fma(x, THIRTYTWO_OVER_PI[2], truncated_prod);
  double k_lo = fputil::nearest_integer(prod_lo);
  y = fputil::fma(x, THIRTYTWO_OVER_PI[2], truncated_prod - k_lo);
  y = fputil::fma(x, THIRTYTWO_OVER_PI[3], y);
  y = fputil::fma(x, THIRTYTWO_OVER_PI[4], y);

  return static_cast<int64_t>(k_lo);
}

} // namespace fma

} // namespace LIBC_NAMESPACE

#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H