1 //===-- Utilities for double precision trigonometric functions ------------===// 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 #include "src/__support/FPUtil/FPBits.h" 10 #include "src/__support/FPUtil/ManipulationFunctions.h" 11 #include "src/__support/FPUtil/UInt.h" 12 #include "src/__support/FPUtil/XFloat.h" 13 14 using FPBits = __llvm_libc::fputil::FPBits<double>; 15 16 namespace __llvm_libc { 17 18 // Implementation is based on the Payne and Hanek range reduction algorithm. 19 // The caller should ensure that x is positive. 20 // Consider: 21 // x/y = x * 1/y = I + F 22 // I is the integral part and F the fractional part of the result of the 23 // division operation. Then M = mod(x, y) = F * y. In order to compute M, we 24 // first compute F. We do it by dropping bits from 1/y which would only 25 // contribute integral results in the operation x * 1/y. This helps us get 26 // accurate values of F even when x is a very large number. 27 // 28 // Internal operations are performed at 192 bits of precision. 29 static double mod_impl(double x, const uint64_t y_bits[3], 30 const uint64_t inv_y_bits[20], int y_exponent, 31 int inv_y_exponent) { 32 FPBits bits(x); 33 int exponent = bits.getExponent(); 34 int bit_drop = (exponent - 52) + inv_y_exponent + 1; 35 bit_drop = bit_drop >= 0 ? bit_drop : 0; 36 int word_drop = bit_drop / 64; 37 bit_drop %= 64; 38 fputil::UInt<256> man4; 39 for (size_t i = 0; i < 4; ++i) { 40 man4[3 - i] = inv_y_bits[word_drop + i]; 41 } 42 man4.shift_left(bit_drop); 43 fputil::UInt<192> man_bits; 44 for (size_t i = 0; i < 3; ++i) 45 man_bits[i] = man4[i + 1]; 46 fputil::XFloat<192> result(inv_y_exponent - word_drop * 64 - bit_drop, 47 man_bits); 48 result.mul(x); 49 result.drop_int(); // |result| now holds fractional part of x/y. 50 51 fputil::UInt<192> y_man; 52 for (size_t i = 0; i < 3; ++i) 53 y_man[i] = y_bits[2 - i]; 54 fputil::XFloat<192> y_192(y_exponent, y_man); 55 return result.mul(y_192); 56 } 57 58 static constexpr int TwoPIExponent = 2; 59 60 // The mantissa bits of 2 * PI. 61 // The most signification bits are in the first uint64_t word 62 // and the least signification bits are in the last word. The 63 // first word includes the implicit '1' bit. 64 static constexpr uint64_t TwoPI[] = {0xc90fdaa22168c234, 0xc4c6628b80dc1cd1, 65 0x29024e088a67cc74}; 66 67 static constexpr int InvTwoPIExponent = -3; 68 69 // The mantissa bits of 1/(2 * PI). 70 // The most signification bits are in the first uint64_t word 71 // and the least signification bits are in the last word. The 72 // first word includes the implicit '1' bit. 73 static constexpr uint64_t InvTwoPI[] = { 74 0xa2f9836e4e441529, 0xfc2757d1f534ddc0, 0xdb6295993c439041, 75 0xfe5163abdebbc561, 0xb7246e3a424dd2e0, 0x6492eea09d1921c, 76 0xfe1deb1cb129a73e, 0xe88235f52ebb4484, 0xe99c7026b45f7e41, 77 0x3991d639835339f4, 0x9c845f8bbdf9283b, 0x1ff897ffde05980f, 78 0xef2f118b5a0a6d1f, 0x6d367ecf27cb09b7, 0x4f463f669e5fea2d, 79 0x7527bac7ebe5f17b, 0x3d0739f78a5292ea, 0x6bfb5fb11f8d5d08, 80 0x56033046fc7b6bab, 0xf0cfbc209af4361e}; 81 82 double mod_2pi(double x) { 83 static constexpr double _2pi = 6.283185307179586; 84 if (x < _2pi) 85 return x; 86 return mod_impl(x, TwoPI, InvTwoPI, TwoPIExponent, InvTwoPIExponent); 87 } 88 89 // Returns mod(x, pi/2) 90 double mod_pi_over_2(double x) { 91 static constexpr double pi_over_2 = 1.5707963267948966; 92 if (x < pi_over_2) 93 return x; 94 return mod_impl(x, TwoPI, InvTwoPI, TwoPIExponent - 2, InvTwoPIExponent + 2); 95 } 96 97 // Returns mod(x, pi/4) 98 double mod_pi_over_4(double x) { 99 static constexpr double pi_over_4 = 0.7853981633974483; 100 if (x < pi_over_4) 101 return x; 102 return mod_impl(x, TwoPI, InvTwoPI, TwoPIExponent - 3, InvTwoPIExponent + 3); 103 } 104 105 } // namespace __llvm_libc 106