1//===-- Implementation header for asinpif16 ---------------------*- 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_ASINPIF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINPIF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "hdr/errno_macros.h"
17#include "hdr/fenv_macros.h"
18#include "src/__support/FPUtil/FEnvImpl.h"
19#include "src/__support/FPUtil/FPBits.h"
20#include "src/__support/FPUtil/PolyEval.h"
21#include "src/__support/FPUtil/cast.h"
22#include "src/__support/FPUtil/multiply_add.h"
23#include "src/__support/FPUtil/sqrt.h"
24#include "src/__support/macros/optimization.h"
25
26namespace LIBC_NAMESPACE_DECL {
27
28namespace math {
29
30LIBC_INLINE float16 asinpif16(float16 x) {
31 using FPBits = fputil::FPBits<float16>;
32
33 FPBits xbits(x);
34 bool is_neg = xbits.is_neg();
35 double x_abs = fputil::cast<double>(x: xbits.abs().get_val());
36
37 auto signed_result = [is_neg](auto r) -> auto { return is_neg ? -r : r; };
38
39 if (LIBC_UNLIKELY(x_abs > 1.0)) {
40 // aspinf16(NaN) = NaN
41 if (xbits.is_nan()) {
42 if (xbits.is_signaling_nan()) {
43 fputil::raise_except_if_required(FE_INVALID);
44 return FPBits::quiet_nan().get_val();
45 }
46 return x;
47 }
48
49 // 1 < |x| <= +/-inf
50 fputil::raise_except_if_required(FE_INVALID);
51 fputil::set_errno_if_required(EDOM);
52
53 return FPBits::quiet_nan().get_val();
54 }
55
56 // the coefficients for the polynomial approximation of asin(x)/pi in the
57 // range [0, 0.5] extracted using python-sympy
58 //
59 // Python code to generate the coefficients:
60 // > from sympy import *
61 // > import math
62 // > x = symbols('x')
63 // > print(series(asin(x)/math.pi, x, 0, 21))
64 //
65 // OUTPUT:
66 //
67 // 0.318309886183791*x + 0.0530516476972984*x**3 + 0.0238732414637843*x**5 +
68 // 0.0142102627760621*x**7 + 0.00967087327815336*x**9 +
69 // 0.00712127941391293*x**11 + 0.00552355646848375*x**13 +
70 // 0.00444514782463692*x**15 + 0.00367705242846804*x**17 +
71 // 0.00310721681820837*x**19 + O(x**21)
72 //
73 // it's very accurate in the range [0, 0.5] and has a maximum error of
74 // 0.0000000000000001 in the range [0, 0.5].
75 constexpr double POLY_COEFFS[] = {
76 0x1.45f306dc9c889p-2, // x^1
77 0x1.b2995e7b7b5fdp-5, // x^3
78 0x1.8723a1d588a36p-6, // x^5
79 0x1.d1a452f20430dp-7, // x^7
80 0x1.3ce52a3a09f61p-7, // x^9
81 0x1.d2b33e303d375p-8, // x^11
82 0x1.69fde663c674fp-8, // x^13
83 0x1.235134885f19bp-8, // x^15
84 };
85 // polynomial evaluation using horner's method
86 // work only for |x| in [0, 0.5]
87 auto asinpi_polyeval = [&](double x) -> double {
88 return x * fputil::polyeval(x: x * x, a0: POLY_COEFFS[0], a: POLY_COEFFS[1],
89 a: POLY_COEFFS[2], a: POLY_COEFFS[3], a: POLY_COEFFS[4],
90 a: POLY_COEFFS[5], a: POLY_COEFFS[6], a: POLY_COEFFS[7]);
91 };
92
93 // if |x| <= 0.5:
94 if (LIBC_UNLIKELY(x_abs <= 0.5)) {
95 // Use polynomial approximation of asin(x)/pi in the range [0, 0.5]
96 double result = asinpi_polyeval(fputil::cast<double>(x));
97 return fputil::cast<float16>(x: result);
98 }
99
100 // If |x| > 0.5, we need to use the range reduction method:
101 // y = asin(x) => x = sin(y)
102 // because: sin(a) = cos(pi/2 - a)
103 // therefore:
104 // x = cos(pi/2 - y)
105 // let z = pi/2 - y,
106 // x = cos(z)
107 // because: cos(2a) = 1 - 2 * sin^2(a), z = 2a, a = z/2
108 // therefore:
109 // cos(z) = 1 - 2 * sin^2(z/2)
110 // sin(z/2) = sqrt((1 - cos(z))/2)
111 // sin(z/2) = sqrt((1 - x)/2)
112 // let u = (1 - x)/2
113 // then:
114 // sin(z/2) = sqrt(u)
115 // z/2 = asin(sqrt(u))
116 // z = 2 * asin(sqrt(u))
117 // pi/2 - y = 2 * asin(sqrt(u))
118 // y = pi/2 - 2 * asin(sqrt(u))
119 // y/pi = 1/2 - 2 * asin(sqrt(u))/pi
120 //
121 // Finally, we can write:
122 // asinpi(x) = 1/2 - 2 * asinpi(sqrt(u))
123 // where u = (1 - x) /2
124 // = 0.5 - 0.5 * x
125 // = multiply_add(-0.5, x, 0.5)
126
127 double u = fputil::multiply_add(x: -0.5, y: x_abs, z: 0.5);
128 double asinpi_sqrt_u = asinpi_polyeval(fputil::sqrt<double>(x: u));
129 double result = fputil::multiply_add(x: -2.0, y: asinpi_sqrt_u, z: 0.5);
130
131 return fputil::cast<float16>(x: signed_result(result));
132}
133
134} // namespace math
135
136} // namespace LIBC_NAMESPACE_DECL
137
138#endif // LIBC_TYPES_HAS_FLOAT16
139
140#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINPIF16_H
141