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
rsqrt(double x)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,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