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
17namespace LIBC_NAMESPACE_DECL {
18
19// PI and PI / 2
20static constexpr double M_MATH_PI = 0x1.921fb54442d18p+1;
21static constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
22
23extern double ATAN_COEFFS[17][9];
24
25// Look-up table for atan(k/16) with k = 0..16.
26static 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.
48LIBC_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
69LIBC_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]);
92static 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
99LIBC_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

source code of libc/src/math/generic/inv_trigf_utils.h