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 #include "math.h"
34 #include "math_private.h"
35
36 #define _CC (0x1p57L + 1)
37 #define _ROOT sqrtl
38
39 volatile static const double tiny = 1.e-300;
40 static const double half = 0.5, one = 1.;
41
42 /* Full precision high and low parts of 1 / pi. */
43 static const long double
44 invpihi = 3.18309886183790671537767526745028737e-01L,
45 invpilo = -1.28821588763206006125693864783127482e-35L;
46
47 /*
48 * Prior to the leading multiplication by x^2, the rational approximation
49 * has an absolute minimax error less than 2.86e-38 over the [0x1p-56,0.5]
50 * domain (or log2(error) = -124.7).
51 */
52 static inline long double
__r(long double xs)53 __r(long double xs)
54 {
55 static const long double
56 R0 = 5.30516476972984452562945877908381207e-02L,
57 R1 = -2.60404681812983888677486288898051068e-01L,
58 R2 = 5.43274368732204127849880302812989126e-01L,
59 R3 = -6.27476039646895838848725721148947475e-01L,
60 R4 = 4.37806223811161460670580937913479055e-01L,
61 R5 = -1.88854683672201011042310131486279147e-01L,
62 R6 = 4.94528645546233985859708428238164776e-02L,
63 R7 = -7.38228286468784921885881430313681146e-03L,
64 R8 = 5.47609957562826794808966842094709768e-04L,
65 R9 = -1.44837427671490633843418618030301533e-05L,
66 R10= 1.60087061374239702655150492329207605e-08L,
67 S1 = -5.35851261206434695164324875722635942e+00L,
68 S2 = 1.23839541267281630422432341920191869e+01L,
69 S3 = -1.61473998442126533731585929134597693e+01L,
70 S4 = 1.30442314991246752613457576698550221e+01L,
71 S5 = -6.74685395302389922545194846871599697e+00L,
72 S6 = 2.22958089278066716953022085768628547e+00L,
73 S7 = -4.55313732093612723575589438727154778e-01L,
74 S8 = 5.33391487207251214709627478847204628e-02L,
75 S9 = -3.08343599417750507732125775680767303e-03L,
76 S10= 6.12260957946655623349049151336307244e-05L;
77
78 long double r, s;
79 r = R5 + (R6 + (R7 + (R8 + (R9 + R10 * xs) * xs) *xs) * xs) * xs;
80 r = R0 + (R1 + (R2 + (R3 + (R4 + r * xs) * xs) * xs) * xs) * xs;
81 s = S5 + (S6 + (S7 + (S8 + (S9 + S10 * xs) * xs) *xs) * xs) * xs;
82 s = 1 + (S1 + (S2 + (S3 + (S4 + r * xs) * xs) * xs) * xs) * xs;
83 return (xs * (r / s));
84 }
85
86 long double
asinpil(long double x)87 asinpil(long double x)
88 {
89 long double ax, hi, lo, xh, xl, y, zh, zl;
90
91 if (isnan(x) || isinf(x))
92 return ((x - x) / (x - x));
93
94 ax = fabsl(x);
95
96 if (ax > 1) /* |x| > 1 */
97 return ((x - x) / (x - x));
98
99
100 if (ax <= 0.5) { /* |x| <= 0.5 */
101 if (ax < 0x1p-57L) { /* |x| < 0x1p-57 */
102 if (ax < 0x1p-16340L) { /* |x| < 0x1p-16340 */
103 if (ax == 0)
104 return (x);
105 /* Scale for near subnormal. */
106 ax *= 0x1p114;
107 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
108 y = (hi + lo) * 0x1p-114;
109 } else {
110 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
111 y = hi + lo;
112 }
113 } else {
114 y = __r(ax * ax);
115 _XADD(invpihi, invpilo, y, 0, xh, xl);
116 _XMUL(ax, 0, xh, xl, hi, lo);
117 y = hi + lo;
118 }
119 } else if (ax < 1) { /* |x| < 1 */
120 y = 1 - ax;
121 x = __r(y / 2);
122 _XADD(invpihi, invpilo, x, 0, xh, xl);
123 _SQRT(2 * y, zh, zl);
124 _XMUL(xh, xl, zh, zl, hi, lo);
125 _XADD(half, 0, -hi, -lo, y, x);
126 } else /* |x| == 1 */
127 y = half;
128
129 return (x < 0 ? -y : y);
130 }
131
132 /*
133 * See src/s_asinpi.c for implementation details.
134 */
135
136 long double
acospil(long double x)137 acospil(long double x)
138 {
139 long double ax, hi, lo, xh, xl, y, zh, zl;
140
141 if (isnan(x) || isinf(x))
142 return ((x - x) / (x - x));
143
144 ax = fabsl(x);
145
146 if (ax > 1) /* |x| > 1 */
147 return ((x - x) / (x - x));
148
149 if (ax <= 0.5) { /* |x| <= 0.5 */
150 if (ax < 0x1p-55L) { /* |x| < 0x1p-55 */
151 y = (ax == 0) ? half : (ax < 0x1p-113 ?
152 half - tiny : half - x * invpihi);
153 } else {
154 y = __r(x * x);
155 _XADD(invpihi, invpilo, y, 0, xh, xl);
156 _XMUL(x, 0, xh, xl, hi, lo);
157 _XADD(half, 0, -hi, -lo, y, ax);
158 }
159 } else if (ax < 1) { /* |x| < 1 */
160 y = 1 - ax;
161 ax = __r(y / 2);
162 _XADD(invpihi, invpilo, ax, 0, xh, xl);
163 _SQRT(2 * y, zh, zl);
164 _XMUL(xh, xl, zh, zl, hi, lo);
165 if (x < 0)
166 _XADD(one, 0, -hi, -lo, y, ax);
167 else
168 y = hi + lo;
169 } else /* |x| == 1 */
170 y = half;
171
172 return (y);
173 }
174