1 /*-
2 * SPDX-License-Identifier: BSD-2-Clause
3 *
4 * Copyright (c) 2026 Steven G. Kargl
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code must retain the above copyright
11 * notice unmodified, this list of conditions, and the following
12 * disclaimer.
13 * 2. Redistributions in binary form must reproduce the above copyright
14 * notice, this list of conditions and the following disclaimer in the
15 * documentation and/or other materials provided with the distribution.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
18 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 /**
30 * asinpi(x) = asin(x) / pi Eq. (1)
31 *
32 * The rational approximation for asinpi(x) has the following form:
33 *
34 * x R(x^2)
35 * asinpi(x) = ---- + x * x^2 * ------------ Eq. (2)
36 * pi 1 + S(x^2)
37 *
38 * with x^2 = x * x. Define r(x^2) = x^2 * [R / (1 + S)], one then has
39 *
40 * asinpi(x) = x * [1 / pi + r(x^2)] Eq. (3)
41 *
42 * For |x| << 1, asinpi(x) = x / pi. That is, the 2nd term in the righthand
43 * side of Eq. (2) can be neglected for |x| < 0x1p{-N/2} where {N/2} is half
44 * the precision (e.g., N = 53, {N/2} = 26). It is noted that for
45 * |x| < 0x1p{emin+m} with m chosen through testing, x / pi is approaching
46 * or is subnormal. To compute the result, x is scaled by 0x1p{N+1}, x / pi
47 * is computed, and finally rescaled by 0x1p{-(N+1)}.
48 *
49 * In the domain, 0x1p{-N/2} <= |x| < 0.5, the approximation becomes
50 *
51 * asinpi(x) = x * (lo + r(x^2) + hi)
52 *
53 * where lo and hi are full-precision low and high parts of 1 / pi.
54 *
55 * In the interval [0.5,1), the following relationship
56 *
57 * asin(x) = pi / 2 - 2 * asin(t)
58 *
59 * with t = [(1 - x) / 2]^{1/2} is used to rewrite Eq. (1). Thus,
60 *
61 * asinpi(x) = 1 / 2 - 2 * t * [1 / pi + r(t^2)] Eq. (4)
62 *
63 * Note, special cases:
64 *
65 * asinpi(+-0) = +-0, exactly.
66 * asinpi(+-1) = +-1/2, exactly.
67 * asinpi(x) = nan for |x| > 1
68 * asinpi(nan) = nan
69 */
70
71 #include <float.h>
72
73 #include "math.h"
74 #include "math_private.h"
75
76 #define _CC (0x1p27 + 1)
77 #define _ROOT sqrt
78
79 volatile static const double tiny = 1.e-300;
80 static const double half = 0.5, one = 1.;
81
82 /* Full precision high and low parts of 1 / pi. */
83 static const double
84 invpihi = 3.1830988618379069e-01,
85 invpilo = -1.9678676675182486e-17;
86
87 /*
88 * R(x^2)
89 * __r(x^2) = x^2 * ------------
90 * 1 + S(x^2)
91 *
92 * Prior to the leading multiplication by x^2, the rational approximation
93 * has an absolute minimax error less than 8.57e-21 over the [0x1p-40,0.5]
94 * domain (or log2(error) = -66.7).
95 */
96 static inline double
__r(double xs)97 __r(double xs)
98 {
99 static const double
100 R0 = 5.3051647697298449e-02,
101 R1 = -1.2219903601836109e-01,
102 R2 = 9.7236612309627199e-02,
103 R3 = -3.0778625727037261e-02,
104 R4 = 3.1527637063244254e-03,
105 R5 = -1.9159514282614908e-05,
106 S1 = -2.7533975629862262e+00,
107 S2 = 2.8040387218379421e+00,
108 S3 = -1.2867553139513013e+00,
109 S4 = 2.5507476275412666e-01,
110 S5 = -1.6150977787265989e-02;
111 double r, s;
112 r = R0 + (R1 + (R2 + (R3 + (R4 + R5 * xs) * xs) * xs) * xs) * xs;
113 s = 1 + (S1 + (S2 + (S3 + (S4 + S5 * xs) * xs) * xs) * xs) * xs;
114 return (xs * (r / s));
115 }
116
117 #include <stdio.h>
118 double
asinpi(double x)119 asinpi(double x)
120 {
121 double ax, hi, lo, xh, xl, y, zh, zl;
122 uint32_t hx, ix, lx;
123
124 EXTRACT_WORDS(hx, lx, x);
125 ix = hx & 0x7fffffff;
126
127 if (ix > 0x3ff00000) /* |x| > 1 */
128 return ((x - x) / (x - x));
129
130 INSERT_WORDS(ax, ix, lx);
131
132 if (ix <= 0x3fe00000) { /* |x| <= 0.5 */
133 if (ix < 0x3e400000) { /* |x| < 0x1p-27 */
134 if (ix < 0x00800000) { /* |x| < 0x1p-1015 */
135 if ((ix | lx) == 0)
136 return (x);
137 /* Scale for near subnormal. */
138 ax *= 0x1p54;
139 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
140 y = (hi + lo) * 0x1p-54;
141 } else {
142 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
143 y = hi + lo;
144 }
145 } else {
146 y = __r(ax * ax);
147 _XADD(invpihi, invpilo, y, 0, xh, xl);
148 _XMUL(ax, 0, xh, xl, hi, lo);
149 y = hi + lo;
150 }
151 } else if (ix < 0x3ff00000) { /* |x| < 1 */
152 y = 1 - ax;
153 x = __r(y / 2);
154 _XADD(invpihi, invpilo, x, 0, xh, xl);
155 _SQRT(2 * y, zh, zl);
156 _XMUL(xh, xl, zh, zl, hi, lo);
157 _XADD(half, 0, -hi, -lo, y, x);
158 } else /* |x| == 1 */
159 y = half;
160
161 return ((hx & 0x80000000) ? -y : y);
162 }
163
164 #if LDBL_MANT_DIG == 53
165 __weak_reference(asinpi, asinpil);
166 #endif
167
168 /*
169 * acospi(x) = acos(x) / pi
170 *
171 * The implementation uses two identities:
172 *
173 * acos(x) = pi / 2 - asin(x) Eq. (5)
174 * acos(-|x|) = pi - acos(|x|) Eq. (6)
175 *
176 * and the definitions for asinpi(x) above. Conversion of Eq. (5)
177 * with the aid of Eq. (3) leads to the form:
178 *
179 * acospi(x) = 1 / 2 - x * [1 / pi + r(x^2)] Eq. (7)
180 *
181 * where 0 <= |x| < 0.5. There are two thresholds. For |x| < 0x1p{-N}
182 * acospi(x) = 1/2 - tiny, which raises FE_INEXACT while preventing spurious
183 * underflow. For |x| < 0x1p{-M}, acospi(x) = 1/2 - x / pi where M is
184 * a sloppy threshold determined from testing. For 0.5 <= |x| < 1, there
185 * are two approximations:
186 *
187 * acospi(x) = 2 * t * [1 / pi + r(t^2)] Eq. (8)
188 *
189 * for 0.5 <= x < 1. When -1 < x <= -0.5, the relevant expression is
190 *
191 * acospi(x) = 1 - 2 * t * [1 / pi + r(t^2)] Eq. (9)
192 *
193 * Note, special cases:
194 *
195 * acospi(+-0) = 1/2, exactly
196 * acospi(1) = 0, exactly
197 * acospi(-1) = 1, exactly
198 * asinpi(x) = nan for |x| > 1
199 * asinpi(nan) = nan
200 */
201
202 double
acospi(double x)203 acospi(double x)
204 {
205 double ax, hi, lo, xh, xl, y, zh, zl;
206 uint32_t hx, ix, lx;
207
208 EXTRACT_WORDS(hx, lx, x);
209 ix = hx & 0x7fffffff;
210
211 if (ix > 0x3ff00000) /* |x| > 1 */
212 return ((x - x) / (x - x));
213
214 if (ix <= 0x3fe00000) { /* |x| <= 0.5 */
215 if (ix < 0x3eb00000) { /* |x| < 0x1p-20 */
216 y = ((ix | lx) == 0) ? half : ((ix < 0x3ca00000) ?
217 half - tiny : half - x * invpihi);
218 } else {
219 y = __r(x * x);
220 _XADD(invpihi, invpilo, y, 0, xh, xl);
221 _XMUL(x, 0, xh, xl, hi, lo);
222 _XADD(half, 0, -hi, -lo, y, ax);
223 }
224 } else if (ix < 0x3ff00000) { /* |x| < 1 */
225 INSERT_WORDS(ax, ix, lx);
226 y = 1 - ax;
227 ax = __r(y / 2);
228 _XADD(invpihi, invpilo, ax, 0, xh, xl);
229 _SQRT(2 * y, zh, zl);
230 _XMUL(xh, xl, zh, zl, hi, lo);
231 if (hx & 0x80000000)
232 _XADD(one, 0, -hi, -lo, y, ax);
233 else
234 y = hi + lo;
235 } else /* |x| == 1 */
236 y = hx & 0x80000000 ? 1 : 0;
237
238 return (y);
239 }
240
241 #if LDBL_MANT_DIG == 53
242 __weak_reference(acospi, acospil);
243 #endif
244