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 * atanpi(x) = atan(x) / pi Eq. (1)
31 *
32 * Note, special cases:
33 *
34 * atanpi(+-0) = +-0, exactly.
35 * atanpi(+-inf) = +-1/2, exactly.
36 * atanpi(nan) = nan
37 *
38 * Reflection symmetry atanpi(-|x|) = - atanpi(|x|) allows the
39 * implementation to be defined for x >= 0.
40 *
41 * A rational approximation for atanpi(x) has the following form:
42 *
43 * x R(x^2)
44 * atanpi(x) = ---- + x * x^2 * ------------ Eq. (2)
45 * pi 1 + S(x^2)
46 *
47 * with x^2 = x * x. Define r(x^2) = x^2 * [R / (1 + S)], one then has
48 *
49 * atanpi(x) = x * [1 / pi + r(x^2)] Eq. (3)
50 *
51 * In addition, for some subdomains of x, the addition formula is used.
52 *
53 * atanpi(x) = atanpi(v) + atanpi[(x - v) / (1 + v * x)] Eq. (4)
54 *
55 * In the interval [0,0x1p{-N/2}) with N the precision of the floating
56 * point type, Eq. (2) can be reduced to
57 *
58 * atanpi(x) = x / pi. Eq. (5)
59 *
60 * However, the division by pi (or more appropriately multiplication]
61 * by the reciprocal) causes issues with |x| < 0x1p{emin+m} with 'm'
62 * determined from testing. The result of Eq. (5) approaches or is a
63 * subnormal. Here, x is scaled by 0x1p{N+1}, Eq. (4) is evaluated, and
64 * then the result is scaled by 0x1p{-(N+1)}.
65 *
66 * In the interval [0xp{-N/2}, 0.5], Eq. (2) is evaluated where the
67 * rational approximation has be found by a minimax procedure.
68 *
69 * In the interval [0.5,0.75), the addition formula gives
70 *
71 * atanpi(x) = atanpi(x0) + atanpi[(x - x0)/(1 + x0 * x)] Eq. (6)
72 *
73 * with x0 = 5/8 chosen at the center of the interval.
74 *
75 * In the interval [0.75,1), the addition formula gives
76 *
77 * atanpi(x) = atanpi(x1) + atanpi[(x - x1)/(1 + x1 * x)] Eq. (7)
78 *
79 * with x1 = 7/8 chosen at the center of the interval.
80 *
81 * In the interval [1,2), the addition formula gives
82 *
83 * atanpi(x) = atanpi(x2) + atanpi[(x - x2)/(1 + x2 * x)] Eq. (8)
84 *
85 * with x2 = 1.5 chosen at the center of the interval.
86 *
87 * Finally, in the interval [2,inf) the identity
88 *
89 * atanpi(x) = 1/2 - atanpi(1 / x) Eq. (9)
90 */
91 #include <float.h>
92
93 #include "math.h"
94 #include "math_private.h"
95
96 #define _CC (0x1p27 + 1)
97 #define _ROOT sqrt
98
99 volatile static const double tiny = 1.e-300;
100 static const double half = 0.5, one = 1., qrtr = 0.25;
101 static const double x0 = 0.625, x1 = 0.875, x2 = 1.5;
102
103 /* Full precision high and low parts. */
104 static const double
105 invpihi = 3.1830988618379069e-01, /* 1/pi */
106 invpilo = -1.9678676675182486e-17, /* 1/pi */
107 a0hi = 1.7780768448935275e-01, /* atanpi(x0) */
108 a0lo = 6.7223942595197191e-18, /* atanpi(x0) */
109 a1hi = 2.2881069536505358e-01, /* atanpi(x1) */
110 a1lo = 8.7193139538130510e-18, /* atanpi(x1) */
111 a2hi = 3.1283295818900120e-01, /* atanpi(x2) */
112 a2lo = -1.4076885713501453e-17; /* atanpi(x2) */
113
114 /*
115 * R(x^2)
116 * __r(x^2) = x^2 * ------------
117 * 1 + S(x^2)
118 *
119 * Prior to the leading multiplication by x^2, the rational approximation
120 * has an absolute minimax error less than 6.24e-19 over the [0x1p-40,0.5]
121 * domain (or log2(error) = -63.8).
122 */
123 static inline double
__r(double xs)124 __r(double xs)
125 {
126 static const double
127 R0 = -1.0610329539459690e-01,
128 R1 = -2.0683077993309035e-01,
129 R2 = -1.3099673469163398e-01,
130 R3 = -2.9655125284635996e-02,
131 R4 = -1.7208096636878276e-03,
132 S1 = 2.5493341763221311e+00,
133 S2 = 2.3356442152763948e+00,
134 S3 = 9.2164104384874268e-01,
135 S4 = 1.4526328350834117e-01,
136 S5 = 6.2132943401189099e-03;
137 double r, s;
138 r = R0 + (R1 + (R2 + (R3 + R4 * xs) * xs) * xs) * xs;
139 s = 1 + (S1 + (S2 + (S3 + (S4 + S5 * xs) * xs) * xs) * xs) * xs;
140 return (xs * (r / s));
141 }
142
143 double
atanpi(double x)144 atanpi(double x)
145 {
146 double ax, hi, lo, xh, xl, y, zh, zl;
147 uint32_t hx, ix, lx;
148
149 EXTRACT_WORDS(hx, lx, x);
150 ix = hx & 0x7fffffff;
151
152 /* x = +-inf, nan */
153 if (ix >= 0x7ff00000) {
154 if (ix > 0x7ff00000)
155 return (x + x);
156 return ((hx & 0x80000000) ? -half : half);
157 }
158
159 INSERT_WORDS(ax, ix, lx);
160
161 if (ix <= 0x3fe00000) { /* |x| <= 0.5 */
162 if (ix < 0x3e400000) { /* |x| < 0x1p-27 */
163 if (ix < 0x00800000) { /* |x| < 0x1p-1015 */
164 if ((ix | lx) == 0)
165 return (x);
166 /* Scale for near subnormal. */
167 ax *= 0x1p54;
168 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
169 y = (hi + lo) * 0x1p-54;
170 } else {
171 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
172 y = hi + lo;
173 }
174 } else {
175 y = __r(ax * ax);
176 _XADD(invpihi, invpilo, y, 0, xh, xl);
177 _XMUL(ax, 0, xh, xl, hi, lo);
178 y = hi + lo;
179 }
180 } else if (ix < 0x3ff00000) { /* |x| < 1 */
181 if (ix < 0x3fe80000) { /* |x| < 0.75 */
182 x = (ax - x0) / (1 + x0 * ax);
183 y = __r(x * x);
184 _XADD(invpihi, invpilo, y, 0, xh, xl);
185 _XMUL(x, 0, xh, xl, hi, lo);
186 _XADD(a0hi, a0lo, hi, lo, y, xl);
187 } else {
188 x = (ax - x1) / (1 + x1 * ax);
189 y = __r(x * x);
190 _XADD(invpihi, invpilo, y, 0, xh, xl);
191 _XMUL(x, 0, xh, xl, hi, lo);
192 _XADD(a1hi, a1lo, hi, lo, y, xl);
193 }
194 } else if (ix < 0x40000000) { /* |x| < 2 */
195 if (ix == 0x3ff00000 && lx == 0)
196 return ((hx & 0x80000000) ? -qrtr : qrtr);
197 x = (ax - x2) / (1 + x2 * ax);
198 y = __r(x * x);
199 _XADD(invpihi, invpilo, y, 0, xh, xl);
200 _XMUL(x, 0, xh, xl, hi, lo);
201 _XADD(a2hi, a2lo, hi, lo, y, xl);
202 } else { /* |x| > 2 */
203 x = 1 / ax;
204 y = __r(x * x);
205 _XADD(invpihi, invpilo, y, 0, xh, xl);
206 _XMUL(x, 0, xh, xl, hi, lo);
207 _XADD(half, 0, -hi, -lo, y, x);
208 }
209
210 return ((hx & 0x80000000) ? -y : y);
211 }
212
213 #if LDBL_MANT_DIG == 53
214 __weak_reference(atanpi, atanpil);
215 #endif
216