1 /*- 2 * Copyright (c) 2026 Steven G. Kargl 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright 9 * notice unmodified, this list of conditions, and the following 10 * disclaimer. 11 * 2. Redistributions in binary form must reproduce the above copyright 12 * notice, this list of conditions and the following disclaimer in the 13 * documentation and/or other materials provided with the distribution. 14 * 15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 16 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 17 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 19 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 20 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 21 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 */ 26 27 /** 28 * Compute the inverse sqrt of x, i.e., rsqrt(x) = 1 / sqrt(x). 29 * 30 * First, filter out special cases: 31 * 32 * 1. rsqrt(+-0) = +-inf, and raise FE_DIVBYZERO exception. 33 * 2. rsqrt(nan) = NaN. 34 * 3. rsqrt(+inf) returns +0. 35 * 2. rsqrt(x<0) = NaN, and raises FE_INVALID. 36 * 37 * If x is a subnormal, scale x into the normal range by x*0x1pN; while 38 * recording the exponent of the scale factor N. Split the possibly 39 * scaled x into f*2^n with f in [0.5,1). Set m=n or m=n-N (subnormal). 40 * If n is odd, then set f = f/2 and increase n to n+1. Thus, f is 41 * in [0.25,1) with n even. 42 * 43 * An initial estimate of y = rqrt[f](x) is 1 / sqrt[f](x). Exhaustive 44 * testing of rsqrtf() gave a max ULP of 1.49; while testing 500M x in 45 * [0,1000] gave a max ULP of 1.24 for rsqrt(). The value of y is then 46 * used with one iteration of Goldschmidt's algorithm: 47 * 48 * z = x * y 49 * h = y / 2 50 * r = 0.5 - h * z 51 * y = h * r + h 52 * 53 * A factor of 2 appears missing in the above, but it is included in the 54 * exponent m. 55 */ 56 #include <fenv.h> 57 #include <float.h> 58 #include "math.h" 59 #include "math_private.h" 60 61 #pragma STDC FENV_ACCESS ON 62 63 #ifdef _CC 64 #undef _CC 65 #endif 66 #define _CC (0x1p27 + 1) 67 68 double 69 rsqrt(double x) 70 { 71 volatile static const double vzero = 0; 72 static const double half = 0.5; 73 int hx, m, rnd; 74 uint32_t lx, ux; 75 double h, ph, pl, rh, rl, y, zh, zl; 76 77 EXTRACT_WORDS(hx, lx, x); 78 ux = (uint32_t)hx & 0x7fffffff; 79 80 /* x = +-0. Raise exception. */ 81 if ((ux | lx) == 0) 82 return (1 / x); 83 84 /* x is NaN. */ 85 if (ux > 0x7ff00000) 86 return (x + x); 87 88 /* x is +-inf. */ 89 if (ux == 0x7ff00000) 90 return (hx & 0x80000000 ? vzero / vzero : 0.); 91 92 /* x < 0. Raise exception. */ 93 if (hx < 0) 94 return (vzero / vzero); 95 96 /* 97 * If x is subnormal, then scale it into the normal range. 98 * Split x into significand and exponent, x = f * 2^m, with 99 * f in [0.5,1) and m a biased exponent. 100 */ 101 m = 0; 102 if (hx < 0x00100000) { /* Subnormal */ 103 x *= 0x1p54; 104 GET_HIGH_WORD(hx, x); 105 m = -54; 106 } 107 m += (hx >> 20) - 1022; 108 hx = (hx & 0x000fffff) | 0x3fe00000; 109 SET_HIGH_WORD(x, hx); 110 111 /* m is odd. Put x into [0.25,0.5) and increase m. */ 112 if (m & 1) { 113 x /= 2; 114 m += 1; 115 } 116 m = -(m >> 1); /* Prepare for 2^(-m/2). */ 117 118 y = 1 / sqrt(x); /* ~52-bit estimate. */ 119 120 h = y / 2; 121 122 /* 123 * For values of x with a representation of 0x1.ffffffffffffepN 124 * with N an odd integer, the computed rsqrt() is not correctly 125 * rounded in round-to-nearest without toggling the rounding mode 126 * to FE_TOWARDZERO. Note, FE_DOWNWARD also works. However, 127 * messing with the rounding mode is expensive, so only do it 128 * when necessary. Example, x = 3.9999999999999991 129 * gives y --> hx = 0x3ff00000, lx = 0x00000001 130 */ 131 EXTRACT_WORDS(hx, lx, y); 132 if ((hx & 0x000fffff) == 0 && lx == 1) { 133 rnd = fegetround(); 134 fesetround(FE_TOWARDZERO); 135 _MUL(x, y, zh, zl); 136 _XMUL(zh, zl, h, 0, ph, pl); 137 fesetround(rnd); 138 } else { 139 _MUL(x, y, zh, zl); 140 _XMUL(zh, zl, h, 0, ph, pl); 141 } 142 143 _XADD(-ph, -pl, half, 0, rh, rl); 144 y = rh * h + h; 145 146 ux = (m + 1024) << 20; 147 INSERT_WORDS(x, ux, 0); 148 return (y *= x); 149 } 150 151 #if LDBL_MANT_DIG == 53 152 __weak_reference(rsqrt, rsqrtl); 153 #endif 154