//===-- 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/except_value_utils.h"
#include "src/__support/FPUtil/nearest_integer.h"

namespace __llvm_libc {

namespace fma {

static constexpr uint32_t FAST_PASS_BOUND = 0x5880'0000U; // 2^50

// Digits of 1/pi, generated by Sollya with:
// > a0 = D(1/pi);
// > a1 = D(1/pi - a0);
// > a2 = D(1/pi - a0 - a1);
// > a3 = D(1/pi - a0 - a1 - a2);
static constexpr double ONE_OVER_PI[5] = {
    0x1.45f306dc9c883p-2, -0x1.6b01ec5417056p-56, -0x1.6447e493ad4cep-110,
    0x1.e21c820ff28b2p-164, -0x1.508510ea79237p-219};

// Return k and y, where
//   k = round(x / pi) and y = (x / pi) - k.
// Assume x is non-negative.
static inline int64_t small_range_reduction(double x, double &y) {
  double kd = fputil::nearest_integer(x * ONE_OVER_PI[0]);
  y = fputil::fma(x, ONE_OVER_PI[0], -kd);
  y = fputil::fma(x, ONE_OVER_PI[1], y);
  return static_cast<int64_t>(kd);
}

// Return k and y, where
//   k = round(x / pi) and y = (x / pi) - k.
static inline int64_t large_range_reduction(double x, int x_exp, double &y) {
  // 2^50 <= |x| < 2^104
  if (x_exp < 103) {
    // - When x < 2^104, the unit bit is contained in the full exact product of
    // x * ONE_OVER_PI[0].
    // - When 2^50 <= |x| < 2^55, the unit bit is contained
    // in the last 8 bits of double(x * ONE_OVER_PI[0]).
    // - When |x| >= 2^55, the LSB of double(x * ONE_OVER_PI[0]) is at least 2.
    fputil::FPBits<double> prod_hi(x * ONE_OVER_PI[0]);
    prod_hi.bits &= (x_exp < 55) ? (~0xffULL) : (~0ULL); // |x| < 2^55
    double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi));
    double truncated_prod = fputil::fma(x, ONE_OVER_PI[0], -k_hi);
    double prod_lo = fputil::fma(x, ONE_OVER_PI[1], truncated_prod);
    double k_lo = fputil::nearest_integer(prod_lo);
    y = fputil::fma(x, ONE_OVER_PI[1], truncated_prod - k_lo);
    y = fputil::fma(x, ONE_OVER_PI[2], y);
    y = fputil::fma(x, ONE_OVER_PI[3], y);

    return static_cast<int64_t>(k_lo);
  }

  // - When x >= 2^104, the full exact product of x * ONE_OVER_PI[0] does not
  // contain the unit bit, so we can ignore it completely.
  // - When 2^104 <= |x| < 2^109, the unit bit is contained
  // in the last 8 bits of double(x * ONE_OVER_PI[1]).
  // - When |x| >= 2^109, the LSB of double(x * ONE_OVER_PI[1]) is at least 2.
  fputil::FPBits<double> prod_hi(x * ONE_OVER_PI[1]);
  prod_hi.bits &= (x_exp < 109) ? (~0xffULL) : (~0ULL); // |x| < 2^55
  double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi));
  double truncated_prod = fputil::fma(x, ONE_OVER_PI[1], -k_hi);
  double prod_lo = fputil::fma(x, ONE_OVER_PI[2], truncated_prod);
  double k_lo = fputil::nearest_integer(prod_lo);
  y = fputil::fma(x, ONE_OVER_PI[2], truncated_prod - k_lo);
  y = fputil::fma(x, ONE_OVER_PI[3], y);
  y = fputil::fma(x, ONE_OVER_PI[4], y);

  return static_cast<int64_t>(k_lo);
}

// Exceptional cases.
static constexpr int N_EXCEPT_SMALL = 9;

static constexpr fputil::ExceptionalValues<float, N_EXCEPT_SMALL> SmallExcepts{
    /* inputs */ {
        0x3fa7832a, // x = 0x1.4f0654p0
        0x40171973, // x = 0x1.2e32e6p1
        0x4096cbe4, // x = 0x1.2d97c8p2
        0x433b7490, // x = 0x1.76e92p7
        0x437ce5f1, // x = 0x1.f9cbe2p7
        0x46199998, // x = 0x1.33333p13
        0x474d246f, // x = 0x1.9a48dep15
        0x4afdece4, // x = 0x1.fbd9c8p22
        0x55cafb2a, // x = 0x1.95f654p44
    },
    /* outputs (RZ, RU offset, RD offset, RN offset) */
    {
        {0x3f7741b5, 1, 0, 1}, // x = 0x1.4f0654p0, sin(x) = 0x1.ee836ap-1 (RZ)
        {0x3f34290f, 1, 0, 1}, // x = 0x1.2e32e6p1, sin(x) = 0x1.68521ep-1 (RZ)
        {0xbf7fffff, 0, 1, 1}, // x = 0x1.2d97c8p2, sin(x) = -0x1.fffffep-1 (RZ)
        {0xbf5cce62, 0, 1, 0}, // x = 0x1.76e92p7, sin(x) = -0x1.b99cc4p-1 (RZ)
        {0x3f7fffff, 1, 0, 1}, // x = 0x1.f9cbe2p7, sin(x) = 0x1.fffffep-1 (RZ)
        {0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ)
        {0x3f7fffff, 1, 0, 1}, // x = 0x1.9a48dep15, sin(x) = 0x1.fffffep-1 (RZ)
        {0xbf7fb6e0, 0, 1, 1}, // x = 0x1.fbd9c8p22, sin(x) = -0x1.ff6dcp-1 (RZ)
        {0xbf7e7a16, 0, 1,
         1}, // x = 0x1.95f654p44, sin(x) = -0x1.fcf42cp-1 (RZ)
    }};

static constexpr int N_EXCEPT_LARGE = 5;

static constexpr fputil::ExceptionalValues<float, N_EXCEPT_LARGE> LargeExcepts{
    /* inputs */ {
        0x5ebcfdde, // x = 0x1.79fbbcp62
        0x5fa6eba7, // x = 0x1.4dd74ep64
        0x6386134e, // x = 0x1.0c269cp72
        0x6a1976f1, // x = 0x1.32ede2p85
        0x727669d4, // x = 0x1.ecd3a8p101
    },
    /* outputs (RZ, RU offset, RD offset, RN offset) */
    {
        {0x3f50622d, 1, 0, 0}, // x = 0x1.79fbbcp62, sin(x) = 0x1.a0c45ap-1 (RZ)
        {0xbe52464a, 0, 1,
         0}, // x = 0x1.4dd74ep64, sin(x) = -0x1.a48c94p-3 (RZ)
        {0x3f7cb2e7, 1, 0, 0}, // x = 0x1.0c269cp72, sin(x) = 0x1.f965cep-1 (RZ)
        {0x3f7fffff, 1, 0, 1}, // x = 0x1.32ede2p85, sin(x) = 0x1.fffffep-1 (RZ)
        {0xbf7a781d, 0, 1,
         0}, // x = 0x1.ecd3a8p101, sin(x) = -0x1.f4f038p-1 (RZ)
    }};

} // namespace fma

} // namespace __llvm_libc

#endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H