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 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 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