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 * src/s_asinpi.c for implemenation details.
31 */
32
33 #ifdef __i386__
34 #include <ieeefp.h>
35 #endif
36 #include <stdint.h>
37
38 #include "fpmath.h"
39 #include "math.h"
40 #include "math_private.h"
41
42 #define _CC (0x1p32L + 1)
43 #define _ROOT sqrtl
44
45 volatile static const double tiny = 1.e-300;
46 static const double half = 0.5, one = 1.;
47
48 /* 1/pi split into the leading and trailing 53 bits. */
49 static const double
50 invpihi = 3.1830988618379069e-01,
51 invpilo = -1.9678676675182486e-17;
52
53 /*
54 * Prior to the leading multiplication by x^2, the rational approximation
55 * has an absolute minimax error less than 1.36e-22 over the [0x1p-32,0.5]
56 * domain (or log2(error) = -72.6).
57 */
58 static inline long double
__r(long double xs)59 __r(long double xs)
60 {
61 static const union IEEEl2bits
62 R0u = LD80C(0xd94caf3dbdb01c38, -5, 5.30516476972984452564e-02L),
63 R1u = LD80C(0x8d346599ebe3212a, -3, -1.37895190734499743665e-01L),
64 R2u = LD80C(0x852c1751d7112655, -3, 1.30051006670116064731e-01L),
65 R3u = LD80C(0xdb55697553742657, -5, -5.35482520546937456640e-02L),
66 R4u = LD80C(0x931ba448fbd6e8c3, -7, 8.97875827280130762861e-03L),
67 R5u = LD80C(0xdbbe8e80762881b0, -12, -4.19129108247665906395e-04L),
68 S1u = LD80C(0xc3272074944c6389, 1, -3.04926310906120631911e+00L),
69 S2u = LD80C(0xe390d58f48c1a96d, 1, 3.55571497910116337692e+00L),
70 S3u = LD80C(0xfccb6664d1fbfc39, 0, -1.97495727465499715865e+00L),
71 S4u = LD80C(0x86f4f140e3be416b, -1, 5.27175024358933283746e-01L),
72 S5u = LD80C(0xf218070935962bdb, -5, -5.91049456446391222176e-02L),
73 S6u = LD80C(0xec656d26215b624b, -10, 1.80355985054588957542e-03L);
74
75 #define R0 (R0u.e)
76 #define R1 (R1u.e)
77 #define R2 (R2u.e)
78 #define R3 (R3u.e)
79 #define R4 (R4u.e)
80 #define R5 (R5u.e)
81 #define S1 (S1u.e)
82 #define S2 (S2u.e)
83 #define S3 (S3u.e)
84 #define S4 (S4u.e)
85 #define S5 (S5u.e)
86 #define S6 (S6u.e)
87
88 long double r, s;
89 r = R0 + (R1 + (R2 + (R3 + (R4 + R5 * xs) * xs) * xs) * xs) * xs;
90 s = 1 + (S1 + (S2 + (S3 + (S4 + (S5 + S6 * xs) * xs) * xs) *
91 xs) * xs) * xs;
92 return (xs * (r / s));
93 }
94
95 #define GREATER(a) (ix == a && lx > 0x8000000000000000ull)
96 #define LESSEQ(a) (ix == a && lx <= 0x8000000000000000ull)
97
98 long double
asinpil(long double x)99 asinpil(long double x)
100 {
101 long double ax, hi, lo, xh, xl, y, zh, zl;
102 uint64_t lx;
103 uint16_t hx, ix;
104
105 EXTRACT_LDBL80_WORDS(hx, lx, x);
106 ix = hx & 0x7fff;
107
108 if (ix >= 0x4000 || GREATER(0x3fff)) /* |x| > 1 */
109 return ((x - x) / (x - x));
110
111 ENTERI();
112
113 INSERT_LDBL80_WORDS(ax, ix, lx);
114
115 if (ix < 0x3ffe || LESSEQ(0x3ffe)) { /* |x| <= 0.5 */
116 if (ix < 0x3fde) { /* |x| < 0x1p-33 */
117 if (ix < 0x002b) { /* |x| < 0x1p-16340 */
118 if ((ix | lx) == 0)
119 RETURNI(x);
120 /* Scale for near subnormal. */
121 ax *= 0x1p65;
122 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
123 y = (hi + lo) * 0x1p-65;
124 } else {
125 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
126 y = hi + lo;
127 }
128 } else {
129 y = __r(ax * ax);
130 _XADD(invpihi, invpilo, y, 0, xh, xl);
131 _XMUL(ax, 0, xh, xl, hi, lo);
132 y = hi + lo;
133 }
134 } else if (ix < 0x3fff) { /* |x| < 1 */
135 y = 1 - ax;
136 x = __r(y / 2);
137 _XADD(invpihi, invpilo, x, 0, xh, xl); /* 1 / pi + r(t^2) */
138 _SQRT(2 * y, zh, zl); /* 2 * t */
139 _XMUL(xh, xl, zh, zl, hi, lo);
140 _XADD(half, 0, -hi, -lo, y, x);
141 } else /* |x| == 1 */
142 y = half;
143
144 RETURNI((hx & 0x8000) ? -y : y);
145 }
146
147 /*
148 * See src/s_asinpi.c for implementation details.
149 */
150
151 long double
acospil(long double x)152 acospil(long double x)
153 {
154 long double ax, hi, lo, xh, xl, y, zh, zl;
155 uint64_t lx;
156 uint16_t hx, ix;
157
158 EXTRACT_LDBL80_WORDS(hx, lx, x);
159 ix = hx & 0x7fff;
160
161 if (ix >= 0x4000 || GREATER(0x3fff)) /* |x| > 1 */
162 return ((x - x) / (x - x));
163
164 ENTERI();
165
166 if (ix < 0x3ffe || LESSEQ(0x3ffe)) { /* |x| <= 0.5 */
167 if (ix < 0x3fe9) { /* |x| < 0x1p-22 */
168 y = ((ix | lx) == 0) ? half : (LESSEQ(0x3fbf) ?
169 half - tiny : half - x * invpihi);
170 } else {
171 y = __r(x * x);
172 _XADD(invpihi, invpilo, y, 0, xh, xl);
173 _XMUL(x, 0, xh, xl, hi, lo);
174 _XADD(half, 0, -hi, -lo, y, ax);
175 }
176 } else if (ix < 0x3fff) { /* |x| < 1 */
177 INSERT_LDBL80_WORDS(ax, ix, lx);
178 y = 1 - ax;
179 ax = __r(y / 2);
180 _XADD(invpihi, invpilo, ax, 0, xh, xl); /* 1 / pi + r(t^2) */
181 _SQRT(2 * y, zh, zl); /* 2 * t */
182 _XMUL(xh, xl, zh, zl, hi, lo);
183 if (hx & 0x8000)
184 _XADD(one, 0, -hi, -lo, y, ax);
185 else
186 y = hi + lo;
187 } else /* |x| == 1 */
188 y = half;
189
190 RETURNI(y);
191 }
192