libc/src/math/generic/range_reduction.h

   1 //===-- Utilities for trigonometric functions -------------------*- 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_H
  10 #define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_H
  11
  12 #include "src/__support/FPUtil/FPBits.h"
  13 #include "src/__support/FPUtil/multiply_add.h"
  14 #include "src/__support/FPUtil/nearest_integer.h"
  15 #include "src/__support/common.h"
  16
  17 namespace LIBC_NAMESPACE {
  18
  19 namespace generic {
  20
  21 static constexpr uint32_t FAST_PASS_BOUND = 0x4a80'0000U; // 2^22
  22
  23 static constexpr int N_ENTRIES = 8;
  24
  25 // We choose to split bits of 32/pi into 28-bit precision pieces, so that the
  26 // product of x * THIRTYTWO_OVER_PI_28[i] is exact.
  27 // These are generated by Sollya with:
  28 // > a1 = D(round(32/pi, 28, RN)); a1;
  29 // > a2 = D(round(32/pi - a1, 28, RN)); a2;
  30 // > a3 = D(round(32/pi - a1 - a2, 28, RN)); a3;
  31 // > a4 = D(round(32/pi - a1 - a2 - a3, 28, RN)); a4;
  32 // ...
  33 static constexpr double THIRTYTWO_OVER_PI_28[N_ENTRIES] = {
  34     0x1.45f306ep+3,   -0x1.b1bbeaep-28,  0x1.3f84ebp-57,    -0x1.7056592p-87,
  35     0x1.c0db62ap-116, -0x1.4cd8778p-145, -0x1.bef806cp-174, 0x1.63abdecp-204};
  36
  37 // Exponents of the least significant bits of the corresponding entries in
  38 // THIRTYTWO_OVER_PI_28.
  39 static constexpr int THIRTYTWO_OVER_PI_28_LSB_EXP[N_ENTRIES] = {
  40     -24, -55, -81, -114, -143, -170, -200, -230};
  41
  42 // Return k and y, where
  43 //   k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
  44 LIBC_INLINE int64_t small_range_reduction(double x, double &y) {
  45   double prod = x * THIRTYTWO_OVER_PI_28[0];
  46   double kd = fputil::nearest_integer(prod);
  47   y = prod - kd;
  48   y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[1], y);
  49   y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[2], y);
  50   return static_cast<int64_t>(kd);
  51 }
  52
  53 // Return k and y, where
  54 //   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
  55 // For large range, there are at most 2 parts of THIRTYTWO_OVER_PI_28
  56 // contributing to the lowest 6 binary digits (k & 63).  If the least
  57 // significant bit of x * the least significant bit of THIRTYTWO_OVER_PI_28[i]
  58 // >= 64, we can completely ignore THIRTYTWO_OVER_PI_28[i].
  59 LIBC_INLINE int64_t large_range_reduction(double x, int x_exp, double &y) {
  60   int idx = 0;
  61   y = 0;
  62   int x_lsb_exp_m4 = x_exp - fputil::FloatProperties<float>::MANTISSA_WIDTH;
  63
  64   // Skipping the first parts of 32/pi such that:
  65   //   LSB of x * LSB of THIRTYTWO_OVER_PI_28[i] >= 32.
  66   while (x_lsb_exp_m4 + THIRTYTWO_OVER_PI_28_LSB_EXP[idx] > 5)
  67     ++idx;
  68
  69   double prod_hi = x * THIRTYTWO_OVER_PI_28[idx];
  70   // Get the integral part of x * THIRTYTWO_OVER_PI_28[idx]
  71   double k_hi = fputil::nearest_integer(prod_hi);
  72   // Get the fractional part of x * THIRTYTWO_OVER_PI_28[idx]
  73   double frac = prod_hi - k_hi;
  74   double prod_lo = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 1], frac);
  75   double k_lo = fputil::nearest_integer(prod_lo);
  76
  77   // Now y is the fractional parts.
  78   y = prod_lo - k_lo;
  79   y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 2], y);
  80   y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 3], y);
  81
  82   return static_cast<int64_t>(k_hi) + static_cast<int64_t>(k_lo);
  83 }
  84
  85 } // namespace generic
  86
  87 } // namespace LIBC_NAMESPACE
  88
  89 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_H