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 implementation details. 31 */ 32 33 #include "math.h" 34 #include "math_private.h" 35 36 #define _CC (0x1p12F + 1) 37 #define _ROOT sqrtf 38 39 volatile static const float tiny = 1.e-30; 40 static const float half = 0.5f, one = 1.f, qrtr = 0.25f; 41 static const float x0 = 0.625, x1 = 0.875, x2 = 1.5; 42 43 /* Full precision high and low parts. */ 44 static const float 45 invpihi = 3.18309873e-01f, /* 1 / pi */ 46 invpilo = 1.28412765e-08f, /* 1 / pi */ 47 a0hi = 1.77807689e-01f, /* atanpi(x0) */ 48 a0lo = -4.22372093e-09f, /* atanpi(x0) */ 49 a1hi = 2.28810698e-01f, /* atanpi(x1) */ 50 a1lo = -2.42890708e-09f, /* atanpi(x1) */ 51 a2hi = 3.12832952e-01f, /* atanpi(x2) */ 52 a2lo = 6.64328592e-09f; /* atanpi(x2) */ 53 54 /* 55 * Prior to the leading multiplication by x^2, the rational approximation 56 * has an absolute minimax error less than 1.59e-10 over the [0x1p-12,0.5] 57 * domain (or log2(error) = -32.5). 58 */ 59 static inline float 60 __r(float xs) 61 { 62 static const float 63 R0 = -1.06103294e-01f, 64 R1 = -6.81197494e-02f, 65 R2 = -1.61480496e-03f, 66 S1 = 1.24201322e+00f, 67 S2 = 3.31868112e-01f; 68 float r, s; 69 r = R0 + (R1 + R2 * xs) * xs; 70 s = 1 + (S1 + S2 * xs) * xs; 71 return (xs * (r / s)); 72 } 73 74 float 75 atanpif(float x) 76 { 77 float ax, hi, lo, xh, xl, y, zh, zl; 78 uint32_t hx, ix; 79 80 GET_FLOAT_WORD(hx, x); 81 ix = hx & 0x7fffffff; 82 83 /* x = +-inf, nan */ 84 if (ix >= 0x7f800000) { 85 if (ix > 0x7f800000) 86 return (x + x); 87 return ((hx & 0x80000000) ? -half : half); 88 } 89 90 SET_FLOAT_WORD(ax, ix); 91 92 if (ix <= 0x3f000000) { /* |x| <= 0.5 */ 93 if (ix < 0x39000000) { /* |x| < 0x1p-13 */ 94 if (ix < 0x03800000) { /* |x| < 0x1p-120 */ 95 if (ix == 0) 96 return (x); 97 /* Scale for near subnormal. */ 98 ax *= 0x1p25f; 99 _XMUL(ax, 0, invpihi, invpilo, hi, lo); 100 y = (hi + lo) * 0x1p-25f; 101 } else { 102 _XMUL(ax, 0, invpihi, invpilo, hi, lo); 103 y = hi + lo; 104 } 105 } else { 106 y = __r(ax * ax); 107 _XADD(invpihi, invpilo, y, 0, xh, xl); 108 _XMUL(ax, 0, xh, xl, hi, lo); 109 y = hi + lo; 110 } 111 } else if (ix < 0x3f800000) { /* |x| < 1 */ 112 if (ix < 0x3f400000) { /* |x| < 0.75 */ 113 x = (ax - x0) / (1 + x0 * ax); 114 y = __r(x * x); 115 _XADD(invpihi, invpilo, y, 0, xh, xl); 116 _XMUL(x, 0, xh, xl, hi, lo); 117 _XADD(a0hi, a0lo, hi, lo, y, xl); 118 } else { /* |x| < 1 */ 119 x = (ax - x1) / (1 + x1 * ax); 120 y = __r(x * x); 121 _XADD(invpihi, invpilo, y, 0, xh, xl); 122 _XMUL(x, 0, xh, xl, hi, lo); 123 _XADD(a1hi, a1lo, hi, lo, y, xl); 124 } 125 } else if (ix < 0x40000000) { /* |x| < 2 */ 126 if (ix == 0x3f800000) 127 return ((hx & 0x80000000) ? -qrtr : qrtr); 128 x = (ax - x2) / (1 + x2 * ax); 129 y = __r(x * x); 130 _XADD(invpihi, invpilo, y, 0, xh, xl); 131 _XMUL(x, 0, xh, xl, hi, lo); 132 _XADD(a2hi, a2lo, hi, lo, y, xl); 133 } else { /* |x| > 2 */ 134 x = 1 / ax; 135 y = __r(x * x); 136 _XADD(invpihi, invpilo, y, 0, xh, xl); 137 _XMUL(x, 0, xh, xl, hi, lo); 138 _XADD(half, 0, -hi, -lo, y, x); 139 } 140 141 return ((hx & 0x80000000) ? -y : y); 142 } 143