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