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/except_value_utils.h" 15 #include "src/__support/FPUtil/nearest_integer.h" 16 17 namespace __llvm_libc { 18 19 namespace fma { 20 21 static constexpr uint32_t FAST_PASS_BOUND = 0x5880'0000U; // 2^50 22 23 // Digits of 1/pi, generated by Sollya with: 24 // > a0 = D(1/pi); 25 // > a1 = D(1/pi - a0); 26 // > a2 = D(1/pi - a0 - a1); 27 // > a3 = D(1/pi - a0 - a1 - a2); 28 static constexpr double ONE_OVER_PI[5] = { 29 0x1.45f306dc9c883p-2, -0x1.6b01ec5417056p-56, -0x1.6447e493ad4cep-110, 30 0x1.e21c820ff28b2p-164, -0x1.508510ea79237p-219}; 31 32 // Return k and y, where 33 // k = round(x / pi) and y = (x / pi) - k. 34 // Assume x is non-negative. 35 static inline int64_t small_range_reduction(double x, double &y) { 36 double kd = fputil::nearest_integer(x * ONE_OVER_PI[0]); 37 y = fputil::fma(x, ONE_OVER_PI[0], -kd); 38 y = fputil::fma(x, ONE_OVER_PI[1], y); 39 return static_cast<int64_t>(kd); 40 } 41 42 // Return k and y, where 43 // k = round(x / pi) and y = (x / pi) - k. 44 static inline int64_t large_range_reduction(double x, int x_exp, double &y) { 45 // 2^50 <= |x| < 2^104 46 if (x_exp < 103) { 47 // - When x < 2^104, the unit bit is contained in the full exact product of 48 // x * ONE_OVER_PI[0]. 49 // - When 2^50 <= |x| < 2^55, the unit bit is contained 50 // in the last 8 bits of double(x * ONE_OVER_PI[0]). 51 // - When |x| >= 2^55, the LSB of double(x * ONE_OVER_PI[0]) is at least 2. 52 fputil::FPBits<double> prod_hi(x * ONE_OVER_PI[0]); 53 prod_hi.bits &= (x_exp < 55) ? (~0xffULL) : (~0ULL); // |x| < 2^55 54 double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi)); 55 double truncated_prod = fputil::fma(x, ONE_OVER_PI[0], -k_hi); 56 double prod_lo = fputil::fma(x, ONE_OVER_PI[1], truncated_prod); 57 double k_lo = fputil::nearest_integer(prod_lo); 58 y = fputil::fma(x, ONE_OVER_PI[1], truncated_prod - k_lo); 59 y = fputil::fma(x, ONE_OVER_PI[2], y); 60 y = fputil::fma(x, ONE_OVER_PI[3], y); 61 62 return static_cast<int64_t>(k_lo); 63 } 64 65 // - When x >= 2^104, the full exact product of x * ONE_OVER_PI[0] does not 66 // contain the unit bit, so we can ignore it completely. 67 // - When 2^104 <= |x| < 2^109, the unit bit is contained 68 // in the last 8 bits of double(x * ONE_OVER_PI[1]). 69 // - When |x| >= 2^109, the LSB of double(x * ONE_OVER_PI[1]) is at least 2. 70 fputil::FPBits<double> prod_hi(x * ONE_OVER_PI[1]); 71 prod_hi.bits &= (x_exp < 109) ? (~0xffULL) : (~0ULL); // |x| < 2^55 72 double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi)); 73 double truncated_prod = fputil::fma(x, ONE_OVER_PI[1], -k_hi); 74 double prod_lo = fputil::fma(x, ONE_OVER_PI[2], truncated_prod); 75 double k_lo = fputil::nearest_integer(prod_lo); 76 y = fputil::fma(x, ONE_OVER_PI[2], truncated_prod - k_lo); 77 y = fputil::fma(x, ONE_OVER_PI[3], y); 78 y = fputil::fma(x, ONE_OVER_PI[4], y); 79 80 return static_cast<int64_t>(k_lo); 81 } 82 83 // Exceptional cases. 84 static constexpr int N_EXCEPT_SMALL = 9; 85 86 static constexpr fputil::ExceptionalValues<float, N_EXCEPT_SMALL> SmallExcepts{ 87 /* inputs */ { 88 0x3fa7832a, // x = 0x1.4f0654p0 89 0x40171973, // x = 0x1.2e32e6p1 90 0x4096cbe4, // x = 0x1.2d97c8p2 91 0x433b7490, // x = 0x1.76e92p7 92 0x437ce5f1, // x = 0x1.f9cbe2p7 93 0x46199998, // x = 0x1.33333p13 94 0x474d246f, // x = 0x1.9a48dep15 95 0x4afdece4, // x = 0x1.fbd9c8p22 96 0x55cafb2a, // x = 0x1.95f654p44 97 }, 98 /* outputs (RZ, RU offset, RD offset, RN offset) */ 99 { 100 {0x3f7741b5, 1, 0, 1}, // x = 0x1.4f0654p0, sin(x) = 0x1.ee836ap-1 (RZ) 101 {0x3f34290f, 1, 0, 1}, // x = 0x1.2e32e6p1, sin(x) = 0x1.68521ep-1 (RZ) 102 {0xbf7fffff, 0, 1, 1}, // x = 0x1.2d97c8p2, sin(x) = -0x1.fffffep-1 (RZ) 103 {0xbf5cce62, 0, 1, 0}, // x = 0x1.76e92p7, sin(x) = -0x1.b99cc4p-1 (RZ) 104 {0x3f7fffff, 1, 0, 1}, // x = 0x1.f9cbe2p7, sin(x) = 0x1.fffffep-1 (RZ) 105 {0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ) 106 {0x3f7fffff, 1, 0, 1}, // x = 0x1.9a48dep15, sin(x) = 0x1.fffffep-1 (RZ) 107 {0xbf7fb6e0, 0, 1, 1}, // x = 0x1.fbd9c8p22, sin(x) = -0x1.ff6dcp-1 (RZ) 108 {0xbf7e7a16, 0, 1, 109 1}, // x = 0x1.95f654p44, sin(x) = -0x1.fcf42cp-1 (RZ) 110 }}; 111 112 static constexpr int N_EXCEPT_LARGE = 5; 113 114 static constexpr fputil::ExceptionalValues<float, N_EXCEPT_LARGE> LargeExcepts{ 115 /* inputs */ { 116 0x5ebcfdde, // x = 0x1.79fbbcp62 117 0x5fa6eba7, // x = 0x1.4dd74ep64 118 0x6386134e, // x = 0x1.0c269cp72 119 0x6a1976f1, // x = 0x1.32ede2p85 120 0x727669d4, // x = 0x1.ecd3a8p101 121 }, 122 /* outputs (RZ, RU offset, RD offset, RN offset) */ 123 { 124 {0x3f50622d, 1, 0, 0}, // x = 0x1.79fbbcp62, sin(x) = 0x1.a0c45ap-1 (RZ) 125 {0xbe52464a, 0, 1, 126 0}, // x = 0x1.4dd74ep64, sin(x) = -0x1.a48c94p-3 (RZ) 127 {0x3f7cb2e7, 1, 0, 0}, // x = 0x1.0c269cp72, sin(x) = 0x1.f965cep-1 (RZ) 128 {0x3f7fffff, 1, 0, 1}, // x = 0x1.32ede2p85, sin(x) = 0x1.fffffep-1 (RZ) 129 {0xbf7a781d, 0, 1, 130 0}, // x = 0x1.ecd3a8p101, sin(x) = -0x1.f4f038p-1 (RZ) 131 }}; 132 133 } // namespace fma 134 135 } // namespace __llvm_libc 136 137 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H 138