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/except_value_utils.h" 14 #include "src/__support/FPUtil/multiply_add.h" 15 #include "src/__support/FPUtil/nearest_integer.h" 16 17 namespace __llvm_libc { 18 19 namespace generic { 20 21 static constexpr uint32_t FAST_PASS_BOUND = 0x4c80'0000U; // 2^26 22 23 static constexpr int N_ENTRIES = 8; 24 25 // We choose to split bits of 1/pi into 28-bit precision pieces, so that the 26 // product of x * ONE_OVER_PI_28[i] is exact. 27 // These are generated by Sollya with: 28 // > a1 = D(round(1/pi, 28, RN)); a1; 29 // > a2 = D(round(1/pi - a1, 28, RN)); a2; 30 // > a3 = D(round(1/pi - a1 - a2, 28, RN)); a3; 31 // > a4 = D(round(1/pi - a1 - a2 - a3, 28, RN)); a4; 32 // ... 33 static constexpr double ONE_OVER_PI_28[N_ENTRIES] = { 34 0x1.45f306ep-2, -0x1.b1bbeaep-33, 0x1.3f84ebp-62, -0x1.7056592p-92, 35 0x1.c0db62ap-121, -0x1.4cd8778p-150, -0x1.bef806cp-179, 0x1.63abdecp-209}; 36 37 // Exponents of the least significant bits of the corresponding entries in 38 // ONE_OVER_PI_28. 39 static constexpr int ONE_OVER_PI_28_LSB_EXP[N_ENTRIES] = { 40 -29, -60, -86, -119, -148, -175, -205, -235}; 41 42 // Return (k mod 2) and y, where 43 // k = round(x / pi) and y = (x / pi) - k. 44 static inline int64_t small_range_reduction(double x, double &y) { 45 double prod = x * ONE_OVER_PI_28[0]; 46 double kd = fputil::nearest_integer(prod); 47 y = prod - kd; 48 y = fputil::multiply_add(x, ONE_OVER_PI_28[1], y); 49 y = fputil::multiply_add(x, ONE_OVER_PI_28[2], y); 50 return static_cast<int64_t>(kd); 51 } 52 53 // Return k and y, where 54 // k = round(x / pi) and y = (x / pi) - k. 55 // For large range, there are at most 2 parts of ONE_OVER_PI_28 contributing to 56 // the unit binary digit (k & 1). If the least significant bit of x * the least 57 // significant bit of ONE_OVER_PI_28[i] > 1, we can completely ignore 58 // ONE_OVER_PI_28[i]. 59 static 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 = x_exp - fputil::FloatProperties<float>::MANTISSA_WIDTH; 63 64 // Skipping the first parts of 1/pi such that: 65 // LSB of x * LSB of ONE_OVER_PI_28[i] > 1. 66 while (x_lsb_exp + ONE_OVER_PI_28_LSB_EXP[idx] > 0) 67 ++idx; 68 69 double prod_hi = x * ONE_OVER_PI_28[idx]; 70 // Get the integral part of x * ONE_OVER_PI_28[idx] 71 double k_hi = fputil::nearest_integer(prod_hi); 72 // Get the fractional part of x * ONE_OVER_PI_28[idx] 73 double frac = prod_hi - k_hi; 74 double prod_lo = fputil::multiply_add(x, ONE_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, ONE_OVER_PI_28[idx + 2], y); 80 y = fputil::multiply_add(x, ONE_OVER_PI_28[idx + 3], y); 81 82 return static_cast<int64_t>(k_hi + k_lo); 83 } 84 85 // Exceptional cases. 86 static constexpr int N_EXCEPT_SMALL = 4; 87 88 static constexpr fputil::ExceptionalValues<float, N_EXCEPT_SMALL> SmallExcepts{ 89 /* inputs */ { 90 0x3fa7832a, // x = 0x1.4f0654p0 91 0x46199998, // x = 0x1.33333p13 92 0x4afdece4, // x = 0x1.fbd9c8p22 93 0x4c2332e9, // x = 0x1.4665d2p25 94 }, 95 /* outputs (RZ, RU offset, RD offset, RN offset) */ 96 { 97 {0x3f7741b5, 1, 0, 1}, // x = 0x1.4f0654p0, sin(x) = 0x1.ee836ap-1 (RZ) 98 {0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ) 99 {0xbf7fb6e0, 0, 1, 1}, // x = 0x1.fbd9c8p22, sin(x) = -0x1.ff6dcp-1 (RZ) 100 {0xbf7fffff, 0, 1, 101 1}, // x = 0x1.4665d2p25, sin(x) = -0x1.fffffep-1 (RZ) 102 }}; 103 104 static constexpr int N_EXCEPT_LARGE = 5; 105 106 static constexpr fputil::ExceptionalValues<float, N_EXCEPT_LARGE> LargeExcepts{ 107 /* inputs */ { 108 0x523947f6, // x = 0x1.728fecp37 109 0x53b146a6, // x = 0x1.628d4cp40 110 0x55cafb2a, // x = 0x1.95f654p44 111 0x6a1976f1, // x = 0x1.32ede2p85 112 0x77584625, // x = 0x1.b08c4ap111 113 }, 114 /* outputs (RZ, RU offset, RD offset, RN offset) */ 115 { 116 {0xbf12791d, 0, 1, 117 1}, // x = 0x1.728fecp37, sin(x) = -0x1.24f23ap-1 (RZ) 118 {0xbf7fffff, 0, 1, 119 1}, // x = 0x1.628d4cp40, sin(x) = -0x1.fffffep-1 (RZ) 120 {0xbf7e7a16, 0, 1, 121 1}, // x = 0x1.95f654p44, sin(x) = -0x1.fcf42cp-1 (RZ) 122 {0x3f7fffff, 1, 0, 1}, // x = 0x1.32ede2p85, sin(x) = 0x1.fffffep-1 (RZ) 123 {0xbf7fffff, 0, 1, 124 1}, // x = 0x1.b08c4ap111, sin(x) = -0x1.fffffep-1 (RZ) 125 }}; 126 127 } // namespace generic 128 129 } // namespace __llvm_libc 130 131 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_H 132