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