1 | //===-- Single-precision general inverse 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 | #ifndef LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H |
10 | #define LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H |
11 | |
12 | #include "src/__support/FPUtil/PolyEval.h" |
13 | #include "src/__support/FPUtil/multiply_add.h" |
14 | #include "src/__support/common.h" |
15 | #include "src/__support/macros/config.h" |
16 | |
17 | namespace LIBC_NAMESPACE_DECL { |
18 | |
19 | // PI and PI / 2 |
20 | static constexpr double M_MATH_PI = 0x1.921fb54442d18p+1; |
21 | static constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0; |
22 | |
23 | extern double ATAN_COEFFS[17][9]; |
24 | |
25 | // Look-up table for atan(k/16) with k = 0..16. |
26 | static constexpr double ATAN_K_OVER_16[17] = { |
27 | 0.0, |
28 | 0x1.ff55bb72cfdeap-5, |
29 | 0x1.fd5ba9aac2f6ep-4, |
30 | 0x1.7b97b4bce5b02p-3, |
31 | 0x1.f5b75f92c80ddp-3, |
32 | 0x1.362773707ebccp-2, |
33 | 0x1.6f61941e4def1p-2, |
34 | 0x1.a64eec3cc23fdp-2, |
35 | 0x1.dac670561bb4fp-2, |
36 | 0x1.0657e94db30dp-1, |
37 | 0x1.1e00babdefeb4p-1, |
38 | 0x1.345f01cce37bbp-1, |
39 | 0x1.4978fa3269ee1p-1, |
40 | 0x1.5d58987169b18p-1, |
41 | 0x1.700a7c5784634p-1, |
42 | 0x1.819d0b7158a4dp-1, |
43 | 0x1.921fb54442d18p-1, |
44 | }; |
45 | |
46 | // For |x| <= 1/32 and 0 <= i <= 16, return Q(x) such that: |
47 | // Q(x) ~ (atan(x + i/16) - atan(i/16)) / x. |
48 | LIBC_INLINE static double atan_eval(double x, unsigned i) { |
49 | double x2 = x * x; |
50 | |
51 | double c0 = fputil::multiply_add(x, ATAN_COEFFS[i][2], ATAN_COEFFS[i][1]); |
52 | double c1 = fputil::multiply_add(x, ATAN_COEFFS[i][4], ATAN_COEFFS[i][3]); |
53 | double c2 = fputil::multiply_add(x, ATAN_COEFFS[i][6], ATAN_COEFFS[i][5]); |
54 | double c3 = fputil::multiply_add(x, ATAN_COEFFS[i][8], ATAN_COEFFS[i][7]); |
55 | |
56 | double x4 = x2 * x2; |
57 | double d1 = fputil::multiply_add(x2, c1, c0); |
58 | double d2 = fputil::multiply_add(x2, c3, c2); |
59 | double p = fputil::multiply_add(x4, d2, d1); |
60 | return p; |
61 | } |
62 | |
63 | // Evaluate atan without big lookup table. |
64 | // atan(n/d) - atan(k/16) = atan((n/d - k/16) / (1 + (n/d) * (k/16))) |
65 | // = atan((n - d * k/16)) / (d + n * k/16)) |
66 | // So we let q = (n - d * k/16) / (d + n * k/16), |
67 | // and approximate with Taylor polynomial: |
68 | // atan(q) ~ q - q^3/3 + q^5/5 - q^7/7 + q^9/9 |
69 | LIBC_INLINE static double atan_eval_no_table(double num, double den, |
70 | double k_over_16) { |
71 | double num_r = fputil::multiply_add(den, -k_over_16, num); |
72 | double den_r = fputil::multiply_add(num, k_over_16, den); |
73 | double q = num_r / den_r; |
74 | |
75 | constexpr double ATAN_TAYLOR[] = { |
76 | -0x1.5555555555555p-2, |
77 | 0x1.999999999999ap-3, |
78 | -0x1.2492492492492p-3, |
79 | 0x1.c71c71c71c71cp-4, |
80 | }; |
81 | double q2 = q * q; |
82 | double q3 = q2 * q; |
83 | double q4 = q2 * q2; |
84 | double c0 = fputil::multiply_add(q2, ATAN_TAYLOR[1], ATAN_TAYLOR[0]); |
85 | double c1 = fputil::multiply_add(q2, ATAN_TAYLOR[3], ATAN_TAYLOR[2]); |
86 | double d = fputil::multiply_add(q4, c1, c0); |
87 | return fputil::multiply_add(q3, d, q); |
88 | } |
89 | |
90 | // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|], |
91 | // [|1, D...|], [0, 0.5]); |
92 | static constexpr double ASIN_COEFFS[10] = { |
93 | 0x1.5555555540fa1p-3, 0x1.333333512edc2p-4, 0x1.6db6cc1541b31p-5, |
94 | 0x1.f1caff324770ep-6, 0x1.6e43899f5f4f4p-6, 0x1.1f847cf652577p-6, |
95 | 0x1.9b60f47f87146p-7, 0x1.259e2634c494fp-6, -0x1.df946fa875ddp-8, |
96 | 0x1.02311ecf99c28p-5}; |
97 | |
98 | // Evaluate P(x^2) - 1, where P(x^2) ~ asin(x)/x |
99 | LIBC_INLINE static double asin_eval(double xsq) { |
100 | double x4 = xsq * xsq; |
101 | double r1 = fputil::polyeval(x4, ASIN_COEFFS[0], ASIN_COEFFS[2], |
102 | ASIN_COEFFS[4], ASIN_COEFFS[6], ASIN_COEFFS[8]); |
103 | double r2 = fputil::polyeval(x4, ASIN_COEFFS[1], ASIN_COEFFS[3], |
104 | ASIN_COEFFS[5], ASIN_COEFFS[7], ASIN_COEFFS[9]); |
105 | return fputil::multiply_add(xsq, r2, r1); |
106 | } |
107 | |
108 | } // namespace LIBC_NAMESPACE_DECL |
109 | |
110 | #endif // LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H |
111 | |