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.
mod_impl(double x,const uint64_t y_bits[3],const uint64_t inv_y_bits[20],int y_exponent,int inv_y_exponent)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.get_exponent();
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
mod_2pi(double x)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)
mod_pi_over_2(double x)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)
mod_pi_over_4(double x)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