1//===-- Implementation header for asinf -------------------------*- 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___SUPPORT_MATH_ASINF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H
11
12#include "inv_trigf_utils.h"
13#include "src/__support/FPUtil/FEnvImpl.h"
14#include "src/__support/FPUtil/FPBits.h"
15#include "src/__support/FPUtil/except_value_utils.h"
16#include "src/__support/FPUtil/multiply_add.h"
17#include "src/__support/FPUtil/sqrt.h"
18#include "src/__support/macros/config.h"
19#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
20#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
21
22namespace LIBC_NAMESPACE_DECL {
23
24namespace math {
25
26LIBC_INLINE float asinf(float x) {
27 using namespace inv_trigf_utils_internal;
28 using FPBits = typename fputil::FPBits<float>;
29
30 FPBits xbits(x);
31 uint32_t x_uint = xbits.uintval();
32 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
33 constexpr double TWO[2] = {-2.0, 2.0};
34 uint32_t x_sign = x_uint >> 31;
35
36 // |x| <= 0.5-ish
37 if (x_abs < 0x3f04'471dU) {
38 // |x| < 0x1.d12edp-12
39 if (LIBC_UNLIKELY(x_abs < 0x39e8'9768U)) {
40 // When |x| < 2^-12, the relative error of the approximation asin(x) ~ x
41 // is:
42 // |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)
43 // = x^2 / 6
44 // < 2^-25
45 // < epsilon(1)/2.
46 // So the correctly rounded values of asin(x) are:
47 // = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,
48 // or (rounding mode = FE_UPWARD and x is
49 // negative),
50 // = x otherwise.
51 // To simplify the rounding decision and make it more efficient, we use
52 // fma(x, 2^-25, x) instead.
53 // An exhaustive test shows that this formula work correctly for all
54 // rounding modes up to |x| < 0x1.d12edp-12.
55 // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
56 // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
57 // |x| < 2^-125. For targets without FMA instructions, we simply use
58 // double for intermediate results as it is more efficient than using an
59 // emulated version of FMA.
60#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
61 return fputil::multiply_add(x, 0x1.0p-25f, x);
62#else
63 double xd = static_cast<double>(x);
64 return static_cast<float>(fputil::multiply_add(x: xd, y: 0x1.0p-25, z: xd));
65#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
66 }
67
68 // For |x| <= 0.5, we approximate asinf(x) by:
69 // asin(x) = x * P(x^2)
70 // Where P(X^2) = Q(X) is a degree-24 minimax even polynomial approximating
71 // asin(x)/x on [0, 0.5] generated by Sollya with:
72 // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20,
73 // 22, 24|], [|1, D...|], [0, 0.5]);
74 // An exhaustive test shows that this approximation works well up to a
75 // little more than 0.5.
76 double xd = static_cast<double>(x);
77 double xsq = xd * xd;
78 double x3 = xd * xsq;
79 double r = asin_eval(xsq);
80 return static_cast<float>(fputil::multiply_add(x: x3, y: r, z: xd));
81 }
82
83 // |x| > 1, return NaNs.
84 if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
85 if (xbits.is_signaling_nan()) {
86 fputil::raise_except_if_required(FE_INVALID);
87 return FPBits::quiet_nan().get_val();
88 }
89
90 if (x_abs <= 0x7f80'0000U) {
91 fputil::set_errno_if_required(EDOM);
92 fputil::raise_except_if_required(FE_INVALID);
93 }
94
95 return FPBits::quiet_nan().get_val();
96 }
97
98 // When |x| > 0.5, we perform range reduction as follow:
99 //
100 // Assume further that 0.5 < x <= 1, and let:
101 // y = asin(x)
102 // We will use the double angle formula:
103 // cos(2y) = 1 - 2 sin^2(y)
104 // and the complement angle identity:
105 // x = sin(y) = cos(pi/2 - y)
106 // = 1 - 2 sin^2 (pi/4 - y/2)
107 // So:
108 // sin(pi/4 - y/2) = sqrt( (1 - x)/2 )
109 // And hence:
110 // pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )
111 // Equivalently:
112 // asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )
113 // Let u = (1 - x)/2, then:
114 // asin(x) = pi/2 - 2 * asin( sqrt(u) )
115 // Moreover, since 0.5 < x <= 1:
116 // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
117 // And hence we can reuse the same polynomial approximation of asin(x) when
118 // |x| <= 0.5:
119 // asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),
120
121 constexpr double M_PI_OVER_4 = -0x1.921fb54442d18p-1;
122
123 xbits.set_sign(Sign::POS);
124 double sign_two = TWO[x_sign]; // sign * (-2)
125 double uf = fputil::multiply_add(x: -0.5f, y: xbits.get_val(), z: 0.5f);
126 double u = static_cast<double>(uf);
127 double c1 = sign_two * fputil::sqrt<double>(x: u);
128 double c2 = fputil::multiply_add(x: sign_two, y: M_PI_OVER_4, z: c1);
129 double c3 = c1 * u;
130
131 double r = asin_eval(xsq: u);
132 return static_cast<float>(fputil::multiply_add(x: c3, y: r, z: c2));
133}
134
135} // namespace math
136
137} // namespace LIBC_NAMESPACE_DECL
138
139#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINF_H
140