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_atanpi.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., qrtr = 0.25;
41 static const double x0 = 0.625, x1 = 0.875, x2 = 1.5;
42
43 /* Full precision high and low parts. */
44 static const long double
45 invpihi = 3.18309886183790671537767526745028737e-01L, /* 1/pi */
46 invpilo = -1.28821588763206006125693864783127482e-35L, /* 1/pi */
47 a0hi = 1.77807684489352753115503587502248118e-01L, /* atanpi(x0) */
48 a0lo = -1.04565241713355884231841538298149212e-35L, /* atanpi(x0) */
49 a1hi = 2.28810695365053587806047702039884546e-01L, /* atanpi(x1) */
50 a1lo = -5.01451892949247960667313017509908719e-37L, /* atanpi(x1) */
51 a2hi = 3.12832958189001183813747252435221446e-01L, /* atanpi(x2) */
52 a2lo = 1.35071436051692259858806390941906106e-36L; /* atanpi(x2) */
53
54 /*
55 * Prior to the leading multiplication by x^2, the rational approximation
56 * has an absolute minimax error less than 2.86e-38 over the [0x1p-56,0.5]
57 * domain (or log2(error) = -124.7).
58 */
59 static inline long double
__r(long double xs)60 __r(long double xs)
61 {
62 static const long double
63 R0 = 5.30516476972984452562945877908381207e-02L,
64 R1 = -2.60404681812983888677486288898051068e-01L,
65 R2 = 5.43274368732204127849880302812989126e-01L,
66 R3 = -6.27476039646895838848725721148947475e-01L,
67 R4 = 4.37806223811161460670580937913479055e-01L,
68 R5 = -1.88854683672201011042310131486279147e-01L,
69 R6 = 4.94528645546233985859708428238164776e-02L,
70 R7 = -7.38228286468784921885881430313681146e-03L,
71 R8 = 5.47609957562826794808966842094709768e-04L,
72 R9 = -1.44837427671490633843418618030301533e-05L,
73 R10= 1.60087061374239702655150492329207605e-08L,
74 S1 = -5.35851261206434695164324875722635942e+00L,
75 S2 = 1.23839541267281630422432341920191869e+01L,
76 S3 = -1.61473998442126533731585929134597693e+01L,
77 S4 = 1.30442314991246752613457576698550221e+01L,
78 S5 = -6.74685395302389922545194846871599697e+00L,
79 S6 = 2.22958089278066716953022085768628547e+00L,
80 S7 = -4.55313732093612723575589438727154778e-01L,
81 S8 = 5.33391487207251214709627478847204628e-02L,
82 S9 = -3.08343599417750507732125775680767303e-03L,
83 S10= 6.12260957946655623349049151336307244e-05L;
84
85 long double r, s;
86 r = R5 + (R6 + (R7 + (R8 + (R9 + R10 * xs) * xs) *xs) * xs) * xs;
87 r = R0 + (R1 + (R2 + (R3 + (R4 + r * xs) * xs) * xs) * xs) * xs;
88 s = S5 + (S6 + (S7 + (S8 + (S9 + S10 * xs) * xs) *xs) * xs) * xs;
89 s = 1 + (S1 + (S2 + (S3 + (S4 + r * xs) * xs) * xs) * xs) * xs;
90 return (xs * (r / s));
91 }
92
93 long double
atanpil(long double x)94 atanpil(long double x)
95 {
96 long double ax, hi, lo, xh, xl, y, zh, zl;
97
98 if (isnan(x) || isinf(x))
99 return ((x - x) / (x - x));
100
101 ax = fabsl(x);
102
103 if (ax > 1) /* |x| > 1 */
104 return ((x - x) / (x - x));
105
106
107 if (ax <= 0.5) { /* |x| <= 0.5 */
108 if (ax < 0x1p-57L) { /* |x| < 0x1p-57 */
109 if (ax < 0x1p-16340L) { /* |x| < 0x1p-16340 */
110 if (ax == 0)
111 return (x);
112 /* Scale for near subnormal. */
113 ax *= 0x1p114;
114 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
115 y = (hi + lo) * 0x1p-114;
116 } else {
117 _XMUL(ax, 0, invpihi, invpilo, hi, lo);
118 y = hi + lo;
119 }
120 } else {
121 y = __r(ax * ax);
122 _XADD(invpihi, invpilo, y, 0, xh, xl);
123 _XMUL(ax, 0, xh, xl, hi, lo);
124 y = hi + lo;
125 }
126 } else if (ax < 1) { /* |x| < 1 */
127 if (ax < 0.75) { /* |x| < 0.75 */
128 x = (ax - x0) / (1 + x0 * ax);
129 y = __r(x * x);
130 _XADD(invpihi, invpilo, y, 0, xh, xl);
131 _XMUL(x, 0, xh, xl, hi, lo);
132 _XADD(a0hi, a0lo, hi, lo, y, xl);
133 } else {
134 x = (ax - x1) / (1 + x1 * ax);
135 y = __r(x * x);
136 _XADD(invpihi, invpilo, y, 0, xh, xl);
137 _XMUL(x, 0, xh, xl, hi, lo);
138 _XADD(a1hi, a1lo, hi, lo, y, xl);
139 }
140 } else if (ax < 2) { /* |x| < 2 */
141 if (ax == 1)
142 return (x < 0 ? -qrtr : qrtr);
143 x = (ax - x2) / (1 + x2 * ax);
144 y = __r(x * x);
145 _XADD(invpihi, invpilo, y, 0, xh, xl);
146 _XMUL(x, 0, xh, xl, hi, lo);
147 _XADD(a2hi, a2lo, hi, lo, y, xl);
148 } else { /* |x| > 2 */
149 x = 1 / ax;
150 y = __r(x * x);
151 _XADD(invpihi, invpilo, y, 0, xh, xl);
152 _XMUL(x, 0, xh, xl, hi, lo);
153 _XADD(half, 0, -hi, -lo, y, x);
154 }
155
156 return (x < 0 ? -y : y);
157 }
158