1
/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
3
* ====================================================
4
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6
* Developed at SunSoft, a Sun Microsystems, Inc. business.
7
* Permission to use, copy, modify, and distribute this
8
* software is freely granted, provided that this notice
10
* ====================================================
13
* Return correctly rounded sqrt.
14
* ------------------------------------------
15
* | Use the hardware sqrt if you have one |
16
* ------------------------------------------
18
* Bit by bit method using integer arithmetic. (Slow, but portable)
20
* Scale x to y in [1,4) with even powers of 2:
21
* find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
22
* sqrt(x) = 2^k * sqrt(y)
23
* 2. Bit by bit computation
24
* Let q = sqrt(y) truncated to i bit after binary point (q = 1),
27
* s = 2*q , and y = 2 * ( y - q ). (1)
30
* To compute q from q , one checks whether
37
* If (2) is false, then q = q ; otherwise q = q + 2 .
40
* With some algebric manipulation, it is not difficult to see
41
* that (2) is equivalent to
46
* The advantage of (3) is that s and y can be computed by
48
* the following recurrence formula:
56
* s = s + 2 , y = y - s - 2 (5)
59
* One may easily use induction to prove (4) and (5).
60
* Note. Since the left hand side of (3) contain only i+2 bits,
61
* it does not necessary to do a full (53-bit) comparison
64
* After generating the 53 bits result, we compute one more bit.
65
* Together with the remainder, we can decide whether the
66
* result is exact, bigger than 1/2ulp, or less than 1/2ulp
67
* (it will never equal to 1/2ulp).
68
* The rounding mode can be detected by checking whether
69
* huge + tiny is equal to huge, and whether huge - tiny is
70
* equal to huge for some floating point number "huge" and "tiny".
73
* sqrt(+-0) = +-0 ... exact
75
* sqrt(-ve) = NaN ... with invalid signal
76
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN
81
static const double tiny = 1.0e-300;
86
int32_t sign = (int)0x80000000;
87
int32_t ix0,s0,q,m,t,i;
88
uint32_t r,t1,s1,ix1,q1;
90
EXTRACT_WORDS(ix0, ix1, x);
92
/* take care of Inf and NaN */
93
if ((ix0&0x7ff00000) == 0x7ff00000) {
94
return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
96
/* take care of zero */
98
if (((ix0&~sign)|ix1) == 0)
99
return x; /* sqrt(+-0) = +-0 */
101
return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
105
if (m == 0) { /* subnormal x */
111
for (i=0; (ix0&0x00100000) == 0; i++)
117
m -= 1023; /* unbias exponent */
118
ix0 = (ix0&0x000fffff)|0x00100000;
119
if (m & 1) { /* odd m, double x to make it even */
120
ix0 += ix0 + ((ix1&sign)>>31);
123
m >>= 1; /* m = [m/2] */
125
/* generate sqrt(x) bit by bit */
126
ix0 += ix0 + ((ix1&sign)>>31);
128
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
129
r = 0x00200000; /* r = moving bit from right to left */
138
ix0 += ix0 + ((ix1&sign)>>31);
147
if (t < ix0 || (t == ix0 && t1 <= ix1)) {
149
if ((t1&sign) == sign && (s1&sign) == 0)
157
ix0 += ix0 + ((ix1&sign)>>31);
162
/* use floating add to find out rounding direction */
163
if ((ix0|ix1) != 0) {
164
z = 1.0 - tiny; /* raise inexact flag */
167
if (q1 == (uint32_t)0xffffffff) {
170
} else if (z > 1.0) {
171
if (q1 == (uint32_t)0xfffffffe)
178
ix0 = (q>>1) + 0x3fe00000;
183
INSERT_WORDS(z, ix0, ix1);