libc/src/math/generic/dp_trig.cpp

//===-- Utilities for double precision trigonometric 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
//
//===----------------------------------------------------------------------===//

#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/ManipulationFunctions.h"
#include "src/__support/FPUtil/UInt.h"
#include "src/__support/FPUtil/XFloat.h"

using FPBits = __llvm_libc::fputil::FPBits<double>;

namespace __llvm_libc {

// Implementation is based on the Payne and Hanek range reduction algorithm.
// The caller should ensure that x is positive.
// Consider:
//   x/y = x * 1/y = I + F
// I is the integral part and F the fractional part of the result of the
// division operation. Then M = mod(x, y) = F * y. In order to compute M, we
// first compute F. We do it by dropping bits from 1/y which would only
// contribute integral results in the operation x * 1/y. This helps us get
// accurate values of F even when x is a very large number.
//
// Internal operations are performed at 192 bits of precision.
static double mod_impl(double x, const uint64_t y_bits[3],
                       const uint64_t inv_y_bits[20], int y_exponent,
                       int inv_y_exponent) {
  FPBits bits(x);
  int exponent = bits.getExponent();
  int bit_drop = (exponent - 52) + inv_y_exponent + 1;
  bit_drop = bit_drop >= 0 ? bit_drop : 0;
  int word_drop = bit_drop / 64;
  bit_drop %= 64;
  fputil::UInt<256> man4;
  for (size_t i = 0; i < 4; ++i) {
    man4[3 - i] = inv_y_bits[word_drop + i];
  }
  man4.shift_left(bit_drop);
  fputil::UInt<192> man_bits;
  for (size_t i = 0; i < 3; ++i)
    man_bits[i] = man4[i + 1];
  fputil::XFloat<192> result(inv_y_exponent - word_drop * 64 - bit_drop,
                             man_bits);
  result.mul(x);
  result.drop_int(); // |result| now holds fractional part of x/y.

  fputil::UInt<192> y_man;
  for (size_t i = 0; i < 3; ++i)
    y_man[i] = y_bits[2 - i];
  fputil::XFloat<192> y_192(y_exponent, y_man);
  return result.mul(y_192);
}

static constexpr int TwoPIExponent = 2;

// The mantissa bits of 2 * PI.
// The most signification bits are in the first uint64_t word
// and the least signification bits are in the last word. The
// first word includes the implicit '1' bit.
static constexpr uint64_t TwoPI[] = {0xc90fdaa22168c234, 0xc4c6628b80dc1cd1,
                                     0x29024e088a67cc74};

static constexpr int InvTwoPIExponent = -3;

// The mantissa bits of 1/(2 * PI).
// The most signification bits are in the first uint64_t word
// and the least signification bits are in the last word. The
// first word includes the implicit '1' bit.
static constexpr uint64_t InvTwoPI[] = {
    0xa2f9836e4e441529, 0xfc2757d1f534ddc0, 0xdb6295993c439041,
    0xfe5163abdebbc561, 0xb7246e3a424dd2e0, 0x6492eea09d1921c,
    0xfe1deb1cb129a73e, 0xe88235f52ebb4484, 0xe99c7026b45f7e41,
    0x3991d639835339f4, 0x9c845f8bbdf9283b, 0x1ff897ffde05980f,
    0xef2f118b5a0a6d1f, 0x6d367ecf27cb09b7, 0x4f463f669e5fea2d,
    0x7527bac7ebe5f17b, 0x3d0739f78a5292ea, 0x6bfb5fb11f8d5d08,
    0x56033046fc7b6bab, 0xf0cfbc209af4361e};

double mod_2pi(double x) {
  static constexpr double _2pi = 6.283185307179586;
  if (x < _2pi)
    return x;
  return mod_impl(x, TwoPI, InvTwoPI, TwoPIExponent, InvTwoPIExponent);
}

// Returns mod(x, pi/2)
double mod_pi_over_2(double x) {
  static constexpr double pi_over_2 = 1.5707963267948966;
  if (x < pi_over_2)
    return x;
  return mod_impl(x, TwoPI, InvTwoPI, TwoPIExponent - 2, InvTwoPIExponent + 2);
}

// Returns mod(x, pi/4)
double mod_pi_over_4(double x) {
  static constexpr double pi_over_4 = 0.7853981633974483;
  if (x < pi_over_4)
    return x;
  return mod_impl(x, TwoPI, InvTwoPI, TwoPIExponent - 3, InvTwoPIExponent + 3);
}

} // namespace __llvm_libc