1 | /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
|
---|
2 | /*
|
---|
3 | * ====================================================
|
---|
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
---|
5 | *
|
---|
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
|
---|
9 | * is preserved.
|
---|
10 | * ====================================================
|
---|
11 | */
|
---|
12 | /* sqrt(x)
|
---|
13 | * Return correctly rounded sqrt.
|
---|
14 | * ------------------------------------------
|
---|
15 | * | Use the hardware sqrt if you have one |
|
---|
16 | * ------------------------------------------
|
---|
17 | * Method:
|
---|
18 | * Bit by bit method using integer arithmetic. (Slow, but portable)
|
---|
19 | * 1. Normalization
|
---|
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),
|
---|
25 | * i 0
|
---|
26 | * i+1 2
|
---|
27 | * s = 2*q , and y = 2 * ( y - q ). (1)
|
---|
28 | * i i i i
|
---|
29 | *
|
---|
30 | * To compute q from q , one checks whether
|
---|
31 | * i+1 i
|
---|
32 | *
|
---|
33 | * -(i+1) 2
|
---|
34 | * (q + 2 ) <= y. (2)
|
---|
35 | * i
|
---|
36 | * -(i+1)
|
---|
37 | * If (2) is false, then q = q ; otherwise q = q + 2 .
|
---|
38 | * i+1 i i+1 i
|
---|
39 | *
|
---|
40 | * With some algebric manipulation, it is not difficult to see
|
---|
41 | * that (2) is equivalent to
|
---|
42 | * -(i+1)
|
---|
43 | * s + 2 <= y (3)
|
---|
44 | * i i
|
---|
45 | *
|
---|
46 | * The advantage of (3) is that s and y can be computed by
|
---|
47 | * i i
|
---|
48 | * the following recurrence formula:
|
---|
49 | * if (3) is false
|
---|
50 | *
|
---|
51 | * s = s , y = y ; (4)
|
---|
52 | * i+1 i i+1 i
|
---|
53 | *
|
---|
54 | * otherwise,
|
---|
55 | * -i -(i+1)
|
---|
56 | * s = s + 2 , y = y - s - 2 (5)
|
---|
57 | * i+1 i i+1 i i
|
---|
58 | *
|
---|
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
|
---|
62 | * in (3).
|
---|
63 | * 3. Final rounding
|
---|
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".
|
---|
71 | *
|
---|
72 | * Special cases:
|
---|
73 | * sqrt(+-0) = +-0 ... exact
|
---|
74 | * sqrt(inf) = inf
|
---|
75 | * sqrt(-ve) = NaN ... with invalid signal
|
---|
76 | * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
|
---|
77 | */
|
---|
78 |
|
---|
79 | #include "libm.h"
|
---|
80 |
|
---|
81 | static const double tiny = 1.0e-300;
|
---|
82 |
|
---|
83 | double sqrt(double x)
|
---|
84 | {
|
---|
85 | double z;
|
---|
86 | int32_t sign = (int)0x80000000;
|
---|
87 | int32_t ix0,s0,q,m,t,i;
|
---|
88 | uint32_t r,t1,s1,ix1,q1;
|
---|
89 |
|
---|
90 | EXTRACT_WORDS(ix0, ix1, x);
|
---|
91 |
|
---|
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 */
|
---|
95 | }
|
---|
96 | /* take care of zero */
|
---|
97 | if (ix0 <= 0) {
|
---|
98 | if (((ix0&~sign)|ix1) == 0)
|
---|
99 | return x; /* sqrt(+-0) = +-0 */
|
---|
100 | if (ix0 < 0)
|
---|
101 | return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
|
---|
102 | }
|
---|
103 | /* normalize x */
|
---|
104 | m = ix0>>20;
|
---|
105 | if (m == 0) { /* subnormal x */
|
---|
106 | while (ix0 == 0) {
|
---|
107 | m -= 21;
|
---|
108 | ix0 |= (ix1>>11);
|
---|
109 | ix1 <<= 21;
|
---|
110 | }
|
---|
111 | for (i=0; (ix0&0x00100000) == 0; i++)
|
---|
112 | ix0<<=1;
|
---|
113 | m -= i - 1;
|
---|
114 | ix0 |= ix1>>(32-i);
|
---|
115 | ix1 <<= i;
|
---|
116 | }
|
---|
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);
|
---|
121 | ix1 += ix1;
|
---|
122 | }
|
---|
123 | m >>= 1; /* m = [m/2] */
|
---|
124 |
|
---|
125 | /* generate sqrt(x) bit by bit */
|
---|
126 | ix0 += ix0 + ((ix1&sign)>>31);
|
---|
127 | ix1 += ix1;
|
---|
128 | q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
|
---|
129 | r = 0x00200000; /* r = moving bit from right to left */
|
---|
130 |
|
---|
131 | while (r != 0) {
|
---|
132 | t = s0 + r;
|
---|
133 | if (t <= ix0) {
|
---|
134 | s0 = t + r;
|
---|
135 | ix0 -= t;
|
---|
136 | q += r;
|
---|
137 | }
|
---|
138 | ix0 += ix0 + ((ix1&sign)>>31);
|
---|
139 | ix1 += ix1;
|
---|
140 | r >>= 1;
|
---|
141 | }
|
---|
142 |
|
---|
143 | r = sign;
|
---|
144 | while (r != 0) {
|
---|
145 | t1 = s1 + r;
|
---|
146 | t = s0;
|
---|
147 | if (t < ix0 || (t == ix0 && t1 <= ix1)) {
|
---|
148 | s1 = t1 + r;
|
---|
149 | if ((t1&sign) == sign && (s1&sign) == 0)
|
---|
150 | s0++;
|
---|
151 | ix0 -= t;
|
---|
152 | if (ix1 < t1)
|
---|
153 | ix0--;
|
---|
154 | ix1 -= t1;
|
---|
155 | q1 += r;
|
---|
156 | }
|
---|
157 | ix0 += ix0 + ((ix1&sign)>>31);
|
---|
158 | ix1 += ix1;
|
---|
159 | r >>= 1;
|
---|
160 | }
|
---|
161 |
|
---|
162 | /* use floating add to find out rounding direction */
|
---|
163 | if ((ix0|ix1) != 0) {
|
---|
164 | z = 1.0 - tiny; /* raise inexact flag */
|
---|
165 | if (z >= 1.0) {
|
---|
166 | z = 1.0 + tiny;
|
---|
167 | if (q1 == (uint32_t)0xffffffff) {
|
---|
168 | q1 = 0;
|
---|
169 | q++;
|
---|
170 | } else if (z > 1.0) {
|
---|
171 | if (q1 == (uint32_t)0xfffffffe)
|
---|
172 | q++;
|
---|
173 | q1 += 2;
|
---|
174 | } else
|
---|
175 | q1 += q1 & 1;
|
---|
176 | }
|
---|
177 | }
|
---|
178 | ix0 = (q>>1) + 0x3fe00000;
|
---|
179 | ix1 = q1>>1;
|
---|
180 | if (q&1)
|
---|
181 | ix1 |= sign;
|
---|
182 | ix0 += m << 20;
|
---|
183 | INSERT_WORDS(z, ix0, ix1);
|
---|
184 | return z;
|
---|
185 | }
|
---|