1 //===-- Single-precision sin function -------------------------------------===//
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/math/sinf.h"
10 #include "src/__support/FPUtil/BasicOperations.h"
11 #include "src/__support/FPUtil/FEnvImpl.h"
12 #include "src/__support/FPUtil/FPBits.h"
13 #include "src/__support/FPUtil/PolyEval.h"
14 #include "src/__support/FPUtil/except_value_utils.h"
15 #include "src/__support/FPUtil/multiply_add.h"
16 #include "src/__support/common.h"
17 
18 #include <errno.h>
19 
20 #if defined(LIBC_TARGET_HAS_FMA)
21 #include "range_reduction_fma.h"
22 // using namespace __llvm_libc::fma;
23 using __llvm_libc::fma::FAST_PASS_BOUND;
24 using __llvm_libc::fma::large_range_reduction;
25 using __llvm_libc::fma::LargeExcepts;
26 using __llvm_libc::fma::N_EXCEPT_LARGE;
27 using __llvm_libc::fma::N_EXCEPT_SMALL;
28 using __llvm_libc::fma::small_range_reduction;
29 using __llvm_libc::fma::SmallExcepts;
30 #else
31 #include "range_reduction.h"
32 // using namespace __llvm_libc::generic;
33 using __llvm_libc::generic::FAST_PASS_BOUND;
34 using __llvm_libc::generic::large_range_reduction;
35 using __llvm_libc::generic::LargeExcepts;
36 using __llvm_libc::generic::N_EXCEPT_LARGE;
37 using __llvm_libc::generic::N_EXCEPT_SMALL;
38 using __llvm_libc::generic::small_range_reduction;
39 using __llvm_libc::generic::SmallExcepts;
40 #endif
41 
42 namespace __llvm_libc {
43 
44 LLVM_LIBC_FUNCTION(float, sinf, (float x)) {
45   using FPBits = typename fputil::FPBits<float>;
46   FPBits xbits(x);
47 
48   uint32_t x_u = xbits.uintval();
49   uint32_t x_abs = x_u & 0x7fff'ffffU;
50   double xd, y;
51 
52   // Range reduction:
53   // For |x| > pi/16, we perform range reduction as follows:
54   // Find k and y such that:
55   //   x = (k + y) * pi
56   //   k is an integer
57   //   |y| < 0.5
58   // For small range (|x| < 2^50 when FMA instructions are available, 2^26
59   // otherwise), this is done by performing:
60   //   k = round(x * 1/pi)
61   //   y = x * 1/pi - k
62   // For large range, we will omit all the higher parts of 1/pi such that the
63   // least significant bits of their full products with x are larger than 1,
64   // since sin(x + i * 2pi) = sin(x).
65   //
66   // When FMA instructions are not available, we store the digits of 1/pi in
67   // chunks of 28-bit precision.  This will make sure that the products:
68   //   x * ONE_OVER_PI_28[i] are all exact.
69   // When FMA instructions are available, we simply store the digits of 1/pi in
70   // chunks of doubles (53-bit of precision).
71   // So when multiplying by the largest values of single precision, the
72   // resulting output should be correct up to 2^(-208 + 128) ~ 2^-80.  By the
73   // worst-case analysis of range reduction, |y| >= 2^-38, so this should give
74   // us more than 40 bits of accuracy. For the worst-case estimation of range
75   // reduction, see for instances:
76   //   Elementary Functions by J-M. Muller, Chapter 11,
77   //   Handbook of Floating-Point Arithmetic by J-M. Muller et. al.,
78   //   Chapter 10.2.
79   //
80   // Once k and y are computed, we then deduce the answer by the sine of sum
81   // formula:
82   //   sin(x) = sin((k + y)*pi)
83   //          = sin(y*pi) * cos(k*pi) + cos(y*pi) * sin(k*pi)
84   //          = (-1)^(k & 1) * sin(y*pi)
85   //          ~ (-1)^(k & 1) * y * P(y^2)
86   // where y*P(y^2) is a degree-15 minimax polynomial generated by Sollya
87   // with: > Q = fpminimax(sin(x*pi)/x, [|0, 2, 4, 6, 8, 10, 12, 14|],
88   // [|D...|], [0, 0.5]);
89 
90   // |x| <= pi/16
91   if (x_abs <= 0x3e49'0fdbU) {
92     xd = static_cast<double>(x);
93 
94     // |x| < 0x1.d12ed2p-12f
95     if (x_abs < 0x39e8'9769U) {
96       if (unlikely(x_abs == 0U)) {
97         // For signed zeros.
98         return x;
99       }
100       // When |x| < 2^-12, the relative error of the approximation sin(x) ~ x
101       // is:
102       //   |sin(x) - x| / |sin(x)| < |x^3| / (6|x|)
103       //                           = x^2 / 6
104       //                           < 2^-25
105       //                           < epsilon(1)/2.
106       // So the correctly rounded values of sin(x) are:
107       //   = x - sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
108       //                        or (rounding mode = FE_UPWARD and x is
109       //                        negative),
110       //   = x otherwise.
111       // To simplify the rounding decision and make it more efficient, we use
112       //   fma(x, -2^-25, x) instead.
113       // An exhaustive test shows that this formula work correctly for all
114       // rounding modes up to |x| < 0x1.c555dep-11f.
115       // Note: to use the formula x - 2^-25*x to decide the correct rounding, we
116       // do need fma(x, -2^-25, x) to prevent underflow caused by -2^-25*x when
117       // |x| < 2^-125. For targets without FMA instructions, we simply use
118       // double for intermediate results as it is more efficient than using an
119       // emulated version of FMA.
120 #if defined(LIBC_TARGET_HAS_FMA)
121       return fputil::multiply_add(x, -0x1.0p-25f, x);
122 #else
123       return static_cast<float>(fputil::multiply_add(xd, -0x1.0p-25, xd));
124 #endif // LIBC_TARGET_HAS_FMA
125     }
126 
127     // |x| < pi/16.
128     double xsq = xd * xd;
129 
130     // Degree-9 polynomial approximation:
131     //   sin(x) ~ x + a_3 x^3 + a_5 x^5 + a_7 x^7 + a_9 x^9
132     //          = x (1 + a_3 x^2 + ... + a_9 x^8)
133     //          = x * P(x^2)
134     // generated by Sollya with the following commands:
135     // > display = hexadecimal;
136     // > Q = fpminimax(sin(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/16]);
137     double result =
138         fputil::polyeval(xsq, 1.0, -0x1.55555555554c6p-3, 0x1.1111111085e65p-7,
139                          -0x1.a019f70fb4d4fp-13, 0x1.718d179815e74p-19);
140     return xd * result;
141   }
142 
143   bool x_sign = xbits.get_sign();
144 
145   int64_t k;
146   xd = static_cast<double>(x);
147 
148   if (x_abs < FAST_PASS_BOUND) {
149     using ExceptChecker =
150         typename fputil::ExceptionChecker<float, N_EXCEPT_SMALL>;
151     {
152       float result;
153       if (ExceptChecker::check_odd_func(SmallExcepts, x_abs, x_sign, result)) {
154         return result;
155       }
156     }
157 
158     k = small_range_reduction(xd, y);
159   } else {
160     // x is inf or nan.
161     if (unlikely(x_abs >= 0x7f80'0000U)) {
162       if (x_abs == 0x7f80'0000U) {
163         errno = EDOM;
164         fputil::set_except(FE_INVALID);
165       }
166       return x +
167              FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
168     }
169 
170     using ExceptChecker =
171         typename fputil::ExceptionChecker<float, N_EXCEPT_LARGE>;
172     {
173       float result;
174       if (ExceptChecker::check_odd_func(LargeExcepts, x_abs, x_sign, result))
175         return result;
176     }
177 
178     k = large_range_reduction(xd, xbits.get_exponent(), y);
179   }
180 
181   // After range reduction, k = round(x / pi) and y = (x/pi) - k.
182   // So k is an integer and -0.5 <= y <= 0.5.
183   // Then sin(x) = sin(y*pi + k*pi)
184   //             = (-1)^(k & 1) * sin(y*pi)
185   //             ~ (-1)^(k & 1) * y * P(y^2)
186   // where y*P(y^2) is a degree-15 minimax polynomial generated by Sollya
187   // with: > P = fpminimax(sin(x*pi)/x, [|0, 2, 4, 6, 8, 10, 12, 14|],
188   // [|D...|], [0, 0.5]);
189 
190   constexpr double SIGN[2] = {1.0, -1.0};
191 
192   double ysq = y * y;
193   double result =
194       y * fputil::polyeval(ysq, 0x1.921fb54442d17p1, -0x1.4abbce625bd4bp2,
195                            0x1.466bc67750a3fp1, -0x1.32d2cce1612b5p-1,
196                            0x1.507832417bce6p-4, -0x1.e3062119b6071p-8,
197                            0x1.e89c7aa14122dp-12, -0x1.625b1709dece6p-16);
198 
199   return SIGN[k & 1] * result;
200   // }
201 }
202 
203 } // namespace __llvm_libc
204