libc/src/math/generic/range_reduction_fma.h

   1 //===-- Utilities for trigonometric functions with FMA ----------*- C++ -*-===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 #ifndef LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H
  10 #define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H
  11
  12 #include "src/__support/FPUtil/FMA.h"
  13 #include "src/__support/FPUtil/FPBits.h"
  14 #include "src/__support/FPUtil/nearest_integer.h"
  15 #include "src/__support/common.h"
  16
  17 namespace __llvm_libc {
  18
  19 namespace fma {
  20
  21 static constexpr uint32_t FAST_PASS_BOUND = 0x5600'0000U; // 2^45
  22
  23 // Digits of 32/pi, generated by Sollya with:
  24 // > a0 = D(32/pi);
  25 // > a1 = D(32/pi - a0);
  26 // > a2 = D(32/pi - a0 - a1);
  27 // > a3 = D(32/pi - a0 - a1 - a2);
  28 static constexpr double THIRTYTWO_OVER_PI[5] = {
  29     0x1.45f306dc9c883p+3, -0x1.6b01ec5417056p-51, -0x1.6447e493ad4cep-105,
  30     0x1.e21c820ff28b2p-159, -0x1.508510ea79237p-214};
  31
  32 // Return k and y, where
  33 //   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
  34 LIBC_INLINE int64_t small_range_reduction(double x, double &y) {
  35   double kd = fputil::nearest_integer(x * THIRTYTWO_OVER_PI[0]);
  36   y = fputil::fma(x, THIRTYTWO_OVER_PI[0], -kd);
  37   y = fputil::fma(x, THIRTYTWO_OVER_PI[1], y);
  38   return static_cast<int64_t>(kd);
  39 }
  40
  41 // Return k and y, where
  42 //   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
  43 // This is used for sinf, cosf, sincosf.
  44 LIBC_INLINE int64_t large_range_reduction(double x, int x_exp, double &y) {
  45   // 2^45 <= |x| < 2^99
  46   if (x_exp < 99) {
  47     // - When x < 2^99, the full exact product of x * THIRTYTWO_OVER_PI[0]
  48     // contains at least one integral bit <= 2^5.
  49     // - When 2^45 <= |x| < 2^55, the lowest 6 unit bits are contained
  50     // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[0]).
  51     // - When |x| >= 2^55, the LSB of double(x * THIRTYTWO_OVER_PI[0]) is at
  52     // least 2^6.
  53     fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[0]);
  54     prod_hi.bits &= (x_exp < 55) ? (~0xfffULL) : (~0ULL); // |x| < 2^55
  55     double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi));
  56     double truncated_prod = fputil::fma(x, THIRTYTWO_OVER_PI[0], -k_hi);
  57     double prod_lo = fputil::fma(x, THIRTYTWO_OVER_PI[1], truncated_prod);
  58     double k_lo = fputil::nearest_integer(prod_lo);
  59     y = fputil::fma(x, THIRTYTWO_OVER_PI[1], truncated_prod - k_lo);
  60     y = fputil::fma(x, THIRTYTWO_OVER_PI[2], y);
  61     y = fputil::fma(x, THIRTYTWO_OVER_PI[3], y);
  62
  63     return static_cast<int64_t>(k_lo);
  64   }
  65
  66   // - When x >= 2^110, the full exact product of x * THIRTYTWO_OVER_PI[0] does
  67   // not contain any of the lowest 6 unit bits, so we can ignore it completely.
  68   // - When 2^99 <= |x| < 2^110, the lowest 6 unit bits are contained
  69   // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[1]).
  70   // - When |x| >= 2^110, the LSB of double(x * THIRTYTWO_OVER_PI[1]) is at
  71   // least 64.
  72   fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[1]);
  73   prod_hi.bits &= (x_exp < 110) ? (~0xfffULL) : (~0ULL); // |x| < 2^110
  74   double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi));
  75   double truncated_prod = fputil::fma(x, THIRTYTWO_OVER_PI[1], -k_hi);
  76   double prod_lo = fputil::fma(x, THIRTYTWO_OVER_PI[2], truncated_prod);
  77   double k_lo = fputil::nearest_integer(prod_lo);
  78   y = fputil::fma(x, THIRTYTWO_OVER_PI[2], truncated_prod - k_lo);
  79   y = fputil::fma(x, THIRTYTWO_OVER_PI[3], y);
  80   y = fputil::fma(x, THIRTYTWO_OVER_PI[4], y);
  81
  82   return static_cast<int64_t>(k_lo);
  83 }
  84
  85 } // namespace fma
  86
  87 } // namespace __llvm_libc
  88
  89 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H