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