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.
small_range_reduction(double x,double & y)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