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/macros/config.h" 16 17 namespace LIBC_NAMESPACE_DECL { 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<double>(x, THIRTYTWO_OVER_PI[0], -kd); 37 y = fputil::fma<double>(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.set_uintval(prod_hi.uintval() & 55 ((x_exp < 55) ? (~0xfffULL) : (~0ULL))); // |x| < 2^55 56 double k_hi = fputil::nearest_integer(prod_hi.get_val()); 57 double truncated_prod = fputil::fma<double>(x, THIRTYTWO_OVER_PI[0], -k_hi); 58 double prod_lo = 59 fputil::fma<double>(x, THIRTYTWO_OVER_PI[1], truncated_prod); 60 double k_lo = fputil::nearest_integer(prod_lo); 61 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[1], truncated_prod - k_lo); 62 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[2], y); 63 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[3], y); 64 65 return static_cast<int64_t>(k_lo); 66 } 67 68 // - When x >= 2^110, the full exact product of x * THIRTYTWO_OVER_PI[0] does 69 // not contain any of the lowest 6 unit bits, so we can ignore it completely. 70 // - When 2^99 <= |x| < 2^110, the lowest 6 unit bits are contained 71 // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[1]). 72 // - When |x| >= 2^110, the LSB of double(x * THIRTYTWO_OVER_PI[1]) is at 73 // least 64. 74 fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[1]); 75 prod_hi.set_uintval(prod_hi.uintval() & 76 ((x_exp < 110) ? (~0xfffULL) : (~0ULL))); // |x| < 2^110 77 double k_hi = fputil::nearest_integer(prod_hi.get_val()); 78 double truncated_prod = fputil::fma<double>(x, THIRTYTWO_OVER_PI[1], -k_hi); 79 double prod_lo = fputil::fma<double>(x, THIRTYTWO_OVER_PI[2], truncated_prod); 80 double k_lo = fputil::nearest_integer(prod_lo); 81 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[2], truncated_prod - k_lo); 82 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[3], y); 83 y = fputil::fma<double>(x, THIRTYTWO_OVER_PI[4], y); 84 85 return static_cast<int64_t>(k_lo); 86 } 87 88 } // namespace fma 89 90 } // namespace LIBC_NAMESPACE_DECL 91 92 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H 93