1//===-- Implementation header for asinpif -----------------------*- 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_ASINPIF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINPIF_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/PolyEval.h"
16#include "src/__support/FPUtil/cast.h"
17#include "src/__support/FPUtil/except_value_utils.h"
18#include "src/__support/FPUtil/multiply_add.h"
19#include "src/__support/FPUtil/sqrt.h"
20#include "src/__support/macros/optimization.h"
21
22namespace LIBC_NAMESPACE_DECL {
23namespace math {
24
25LIBC_INLINE float asinpif(float x) {
26 using FPBits = fputil::FPBits<float>;
27
28 FPBits xbits(x);
29 bool is_neg = xbits.is_neg();
30 double x_abs = fputil::cast<double>(x: xbits.abs().get_val());
31
32 auto signed_result = [is_neg](auto r) -> auto { return is_neg ? -r : r; };
33
34 if (LIBC_UNLIKELY(x_abs > 1.0)) {
35 if (xbits.is_nan()) {
36 if (xbits.is_signaling_nan()) {
37 fputil::raise_except_if_required(FE_INVALID);
38 return FPBits::quiet_nan().get_val();
39 }
40 return x;
41 }
42
43 fputil::raise_except_if_required(FE_INVALID);
44 fputil::set_errno_if_required(EDOM);
45 return FPBits::quiet_nan().get_val();
46 }
47
48 // if |x| <= 0.5:
49 // asinpi(x) = x * (c0 + x^2 * P1(x^2))
50 if (LIBC_UNLIKELY(x_abs <= 0.5)) {
51 double x_d = fputil::cast<double>(x);
52 double v2 = x_d * x_d;
53 double result = x_d * fputil::multiply_add(
54 x: v2, y: inv_trigf_utils_internal::asinpi_eval(v2),
55 z: inv_trigf_utils_internal::ASINPI_COEFFS[0]);
56 return fputil::cast<float>(x: result);
57 }
58
59 // If |x| > 0.5:
60 // asinpi(x) = 0.5 - 2 * sqrt(u) * P(u)
61 // = 0.5 - 2 * sqrt(u) * (c0 + u * P1(u))
62 // = (0.5 - 2*sqrt(u)*ONE_OVER_PI_HI)
63 // - 2*sqrt(u) * (ONE_OVER_PI_LO + DELTA_C0 + u * P1(u))
64 //
65 // where u = (1 - |x|) / 2, and
66 // ONE_OVER_PI_HI + ONE_OVER_PI_LO = 1/pi to ~106 bits
67 // DELTA_C0 = c0 - ONE_OVER_PI_HI
68 //
69 // ONE_OVER_PI_LO + DELTA_C0 is a single precomputed constant:
70 // = ONE_OVER_PI_LO + (c0 - ONE_OVER_PI_HI)
71 // = c0 - (ONE_OVER_PI_HI - ONE_OVER_PI_LO)
72 // = c0 - 1/pi (to ~106 bits)
73 constexpr double ONE_OVER_PI_HI = 0x1.45f306dc9c883p-2;
74 constexpr double ONE_OVER_PI_LO = -0x1.6b01ec5417056p-56;
75 // C0_MINUS_1OVERPI = c0 - 1/pi = DELTA_C0 + ONE_OVER_PI_LO
76 constexpr double C0_MINUS_1OVERPI =
77 (inv_trigf_utils_internal::ASINPI_COEFFS[0] - ONE_OVER_PI_HI) +
78 ONE_OVER_PI_LO;
79
80 double u = fputil::multiply_add(x: -0.5, y: x_abs, z: 0.5);
81 double sqrt_u = fputil::sqrt<double>(x: u);
82 double neg2_sqrt_u = -2.0 * sqrt_u;
83
84 // tail = (c0 - 1/pi) + u * P1(u)
85 double tail = fputil::multiply_add(
86 x: u, y: inv_trigf_utils_internal::asinpi_eval(v2: u), z: C0_MINUS_1OVERPI);
87
88 double result_hi = fputil::multiply_add(x: neg2_sqrt_u, y: ONE_OVER_PI_HI, z: 0.5);
89 double result = fputil::multiply_add(x: tail, y: neg2_sqrt_u, z: result_hi);
90
91 return fputil::cast<float>(x: signed_result(result));
92}
93
94} // namespace math
95} // namespace LIBC_NAMESPACE_DECL
96
97#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINPIF_H
98