1 
/* * Mode: C; tabwidth: 8; indenttabsmode: nil; cbasicoffset: 4 * 
2 
* 
3 
* ***** BEGIN LICENSE BLOCK ***** 
4 
* Version: MPL 1.1/GPL 2.0/LGPL 2.1 
5 
* 
6 
* The contents of this file are subject to the Mozilla Public License Version 
7 
* 1.1 (the "License"); you may not use this file except in compliance with 
8 
* the License. You may obtain a copy of the License at 
9 
* http://www.mozilla.org/MPL/ 
10 
* 
11 
* Software distributed under the License is distributed on an "AS IS" basis, 
12 
* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 
13 
* for the specific language governing rights and limitations under the 
14 
* License. 
15 
* 
16 
* The Original Code is Mozilla Communicator client code, released 
17 
* March 31, 1998. 
18 
* 
19 
* The Initial Developer of the Original Code is 
20 
* Sun Microsystems, Inc. 
21 
* Portions created by the Initial Developer are Copyright (C) 1998 
22 
* the Initial Developer. All Rights Reserved. 
23 
* 
24 
* Contributor(s): 
25 
* 
26 
* Alternatively, the contents of this file may be used under the terms of 
27 
* either of the GNU General Public License Version 2 or later (the "GPL"), 
28 
* or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 
29 
* in which case the provisions of the GPL or the LGPL are applicable instead 
30 
* of those above. If you wish to allow use of your version of this file only 
31 
* under the terms of either the GPL or the LGPL, and not to allow others to 
32 
* use your version of this file under the terms of the MPL, indicate your 
33 
* decision by deleting the provisions above and replace them with the notice 
34 
* and other provisions required by the GPL or the LGPL. If you do not delete 
35 
* the provisions above, a recipient may use your version of this file under 
36 
* the terms of any one of the MPL, the GPL or the LGPL. 
37 
* 
38 
* ***** END LICENSE BLOCK ***** */ 
39 
/* @(#)e_sqrt.c 1.3 95/01/18 */ 
40 
/* 
41 
* ==================================================== 
42 
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 
43 
* 
44 
* Developed at SunSoft, a Sun Microsystems, Inc. business. 
45 
* Permission to use, copy, modify, and distribute this 
46 
* software is freely granted, provided that this notice 
47 
* is preserved. 
48 
* ==================================================== 
49 
*/ 
50 

51 
/* __ieee754_sqrt(x) 
52 
* Return correctly rounded sqrt. 
53 
*  
54 
*  Use the hardware sqrt if you have one  
55 
*  
56 
* Method: 
57 
* Bit by bit method using integer arithmetic. (Slow, but portable) 
58 
* 1. Normalization 
59 
* Scale x to y in [1,4) with even powers of 2: 
60 
* find an integer k such that 1 <= (y=x*2^(2k)) < 4, then 
61 
* sqrt(y) = 2^k * sqrt(x) 
62 
* 2. Bit by bit computation 
63 
* Let q = sqrt(y) truncated to i bit after binary point (q = 1), 
64 
* i 0 
65 
* i+1 2 
66 
* s = 2*q , and y = 2 * ( y  q ). (1) 
67 
* i i i i 
68 
* 
69 
* To compute q from q , one checks whether 
70 
* i+1 i 
71 
* 
72 
* (i+1) 2 
73 
* (q + 2 ) <= y. (2) 
74 
* i 
75 
* (i+1) 
76 
* If (2) is false, then q = q ; otherwise q = q + 2 . 
77 
* i+1 i i+1 i 
78 
* 
79 
* With some algebric manipulation, it is not difficult to see 
80 
* that (2) is equivalent to 
81 
* (i+1) 
82 
* s + 2 <= y (3) 
83 
* i i 
84 
* 
85 
* The advantage of (3) is that s and y can be computed by 
86 
* i i 
87 
* the following recurrence formula: 
88 
* if (3) is false 
89 
* 
90 
* s = s , y = y ; (4) 
91 
* i+1 i i+1 i 
92 
* 
93 
* otherwise, 
94 
* i (i+1) 
95 
* s = s + 2 , y = y  s  2 (5) 
96 
* i+1 i i+1 i i 
97 
* 
98 
* One may easily use induction to prove (4) and (5). 
99 
* Note. Since the left hand side of (3) contain only i+2 bits, 
100 
* it does not necessary to do a full (53bit) comparison 
101 
* in (3). 
102 
* 3. Final rounding 
103 
* After generating the 53 bits result, we compute one more bit. 
104 
* Together with the remainder, we can decide whether the 
105 
* result is exact, bigger than 1/2ulp, or less than 1/2ulp 
106 
* (it will never equal to 1/2ulp). 
107 
* The rounding mode can be detected by checking whether 
108 
* huge + tiny is equal to huge, and whether huge  tiny is 
109 
* equal to huge for some floating point number "huge" and "tiny". 
110 
* 
111 
* Special cases: 
112 
* sqrt(+0) = +0 ... exact 
113 
* sqrt(inf) = inf 
114 
* sqrt(ve) = NaN ... with invalid signal 
115 
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN 
116 
* 
117 
* Other methods : see the appended file at the end of the program below. 
118 
* 
119 
*/ 
120 

121 
#include "fdlibm.h" 
122 

123 
#if defined(_MSC_VER) 
124 
/* Microsoft Compiler */ 
125 
#pragma warning( disable : 4723 ) /* disables potential divide by 0 warning */ 
126 
#endif 
127 

128 
#ifdef __STDC__ 
129 
static const double one = 1.0, tiny=1.0e300; 
130 
#else 
131 
static double one = 1.0, tiny=1.0e300; 
132 
#endif 
133 

134 
#ifdef __STDC__ 
135 
double __ieee754_sqrt(double x) 
136 
#else 
137 
double __ieee754_sqrt(x) 
138 
double x; 
139 
#endif 
140 
{ 
141 
fd_twoints u; 
142 
double z; 
143 
int sign = (int)0x80000000; 
144 
unsigned r,t1,s1,ix1,q1; 
145 
int ix0,s0,q,m,t,i; 
146 

147 
u.d = x; 
148 
ix0 = __HI(u); /* high word of x */ 
149 
ix1 = __LO(u); /* low word of x */ 
150 

151 
/* take care of Inf and NaN */ 
152 
if((ix0&0x7ff00000)==0x7ff00000) { 
153 
return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf 
154 
sqrt(inf)=sNaN */ 
155 
} 
156 
/* take care of zero */ 
157 
if(ix0<=0) { 
158 
if(((ix0&(~sign))ix1)==0) return x;/* sqrt(+0) = +0 */ 
159 
else if(ix0<0) 
160 
return (xx)/(xx); /* sqrt(ve) = sNaN */ 
161 
} 
162 
/* normalize x */ 
163 
m = (ix0>>20); 
164 
if(m==0) { /* subnormal x */ 
165 
while(ix0==0) { 
166 
m = 21; 
167 
ix0 = (ix1>>11); ix1 <<= 21; 
168 
} 
169 
for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1; 
170 
m = i1; 
171 
ix0 = (ix1>>(32i)); 
172 
ix1 <<= i; 
173 
} 
174 
m = 1023; /* unbias exponent */ 
175 
ix0 = (ix0&0x000fffff)0x00100000; 
176 
if(m&1){ /* odd m, double x to make it even */ 
177 
ix0 += ix0 + ((ix1&sign)>>31); 
178 
ix1 += ix1; 
179 
} 
180 
m >>= 1; /* m = [m/2] */ 
181 

182 
/* generate sqrt(x) bit by bit */ 
183 
ix0 += ix0 + ((ix1&sign)>>31); 
184 
ix1 += ix1; 
185 
q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */ 
186 
r = 0x00200000; /* r = moving bit from right to left */ 
187 

188 
while(r!=0) { 
189 
t = s0+r; 
190 
if(t<=ix0) { 
191 
s0 = t+r; 
192 
ix0 = t; 
193 
q += r; 
194 
} 
195 
ix0 += ix0 + ((ix1&sign)>>31); 
196 
ix1 += ix1; 
197 
r>>=1; 
198 
} 
199 

200 
r = sign; 
201 
while(r!=0) { 
202 
t1 = s1+r; 
203 
t = s0; 
204 
if((t<ix0)((t==ix0)&&(t1<=ix1))) { 
205 
s1 = t1+r; 
206 
if(((int)(t1&sign)==sign)&&(s1&sign)==0) s0 += 1; 
207 
ix0 = t; 
208 
if (ix1 < t1) ix0 = 1; 
209 
ix1 = t1; 
210 
q1 += r; 
211 
} 
212 
ix0 += ix0 + ((ix1&sign)>>31); 
213 
ix1 += ix1; 
214 
r>>=1; 
215 
} 
216 

217 
/* use floating add to find out rounding direction */ 
218 
if((ix0ix1)!=0) { 
219 
z = onetiny; /* trigger inexact flag */ 
220 
if (z>=one) { 
221 
z = one+tiny; 
222 
if (q1==(unsigned)0xffffffff) { q1=0; q += 1;} 
223 
else if (z>one) { 
224 
if (q1==(unsigned)0xfffffffe) q+=1; 
225 
q1+=2; 
226 
} else 
227 
q1 += (q1&1); 
228 
} 
229 
} 
230 
ix0 = (q>>1)+0x3fe00000; 
231 
ix1 = q1>>1; 
232 
if ((q&1)==1) ix1 = sign; 
233 
ix0 += (m <<20); 
234 
u.d = z; 
235 
__HI(u) = ix0; 
236 
__LO(u) = ix1; 
237 
z = u.d; 
238 
return z; 
239 
} 
240 

241 
/* 
242 
Other methods (use floatingpoint arithmetic) 
243 
 
244 
(This is a copy of a drafted paper by Prof W. Kahan 
245 
and K.C. Ng, written in May, 1986) 
246 

247 
Two algorithms are given here to implement sqrt(x) 
248 
(IEEE double precision arithmetic) in software. 
249 
Both supply sqrt(x) correctly rounded. The first algorithm (in 
250 
Section A) uses newton iterations and involves four divisions. 
251 
The second one uses reciproot iterations to avoid division, but 
252 
requires more multiplications. Both algorithms need the ability 
253 
to chop results of arithmetic operations instead of round them, 
254 
and the INEXACT flag to indicate when an arithmetic operation 
255 
is executed exactly with no roundoff error, all part of the 
256 
standard (IEEE 7541985). The ability to perform shift, add, 
257 
subtract and logical AND operations upon 32bit words is needed 
258 
too, though not part of the standard. 
259 

260 
A. sqrt(x) by Newton Iteration 
261 

262 
(1) Initial approximation 
263 

264 
Let x0 and x1 be the leading and the trailing 32bit words of 
265 
a floating point number x (in IEEE double format) respectively 
266 

267 
1 11 52 ...widths 
268 
 
269 
x: s e  f  
270 
 
271 
msb lsb msb lsb ...order 
272 

273 

274 
  
275 
x0: s e  f1  x1:  f2  
276 
  
277 

278 
By performing shifts and subtracts on x0 and x1 (both regarded 
279 
as integers), we obtain an 8bit approximation of sqrt(x) as 
280 
follows. 
281 

282 
k := (x0>>1) + 0x1ff80000; 
283 
y0 := k  T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits 
284 
Here k is a 32bit integer and T1[] is an integer array containing 
285 
correction terms. Now magically the floating value of y (y's 
286 
leading 32bit word is y0, the value of its trailing word is 0) 
287 
approximates sqrt(x) to almost 8bit. 
288 

289 
Value of T1: 
290 
static int T1[32]= { 
291 
0, 1024, 3062, 5746, 9193, 13348, 18162, 23592, 
292 
29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215, 
293 
83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581, 
294 
16499, 12183, 8588, 5674, 3403, 1742, 661, 130,}; 
295 

296 
(2) Iterative refinement 
297 

298 
Apply Heron's rule three times to y, we have y approximates 
299 
sqrt(x) to within 1 ulp (Unit in the Last Place): 
300 

301 
y := (y+x/y)/2 ... almost 17 sig. bits 
302 
y := (y+x/y)/2 ... almost 35 sig. bits 
303 
y := y(yx/y)/2 ... within 1 ulp 
304 

305 

306 
Remark 1. 
307 
Another way to improve y to within 1 ulp is: 
308 

309 
y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x) 
310 
y := y  0x00100006 ... almost 18 sig. bits to sqrt(x) 
311 

312 
2 
313 
(xy )*y 
314 
y := y + 2*  ...within 1 ulp 
315 
2 
316 
3y + x 
317 

318 

319 
This formula has one division fewer than the one above; however, 
320 
it requires more multiplications and additions. Also x must be 
321 
scaled in advance to avoid spurious overflow in evaluating the 
322 
expression 3y*y+x. Hence it is not recommended uless division 
323 
is slow. If division is very slow, then one should use the 
324 
reciproot algorithm given in section B. 
325 

326 
(3) Final adjustment 
327 

328 
By twiddling y's last bit it is possible to force y to be 
329 
correctly rounded according to the prevailing rounding mode 
330 
as follows. Let r and i be copies of the rounding mode and 
331 
inexact flag before entering the square root program. Also we 
332 
use the expression y+ulp for the next representable floating 
333 
numbers (up and down) of y. Note that y+ulp = either fixed 
334 
point y+1, or multiply y by nextafter(1,+inf) in chopped 
335 
mode. 
336 

337 
I := FALSE; ... reset INEXACT flag I 
338 
R := RZ; ... set rounding mode to roundtowardzero 
339 
z := x/y; ... chopped quotient, possibly inexact 
340 
If(not I) then { ... if the quotient is exact 
341 
if(z=y) { 
342 
I := i; ... restore inexact flag 
343 
R := r; ... restore rounded mode 
344 
return sqrt(x):=y. 
345 
} else { 
346 
z := z  ulp; ... special rounding 
347 
} 
348 
} 
349 
i := TRUE; ... sqrt(x) is inexact 
350 
If (r=RN) then z=z+ulp ... roundedtonearest 
351 
If (r=RP) then { ... roundtoward+inf 
352 
y = y+ulp; z=z+ulp; 
353 
} 
354 
y := y+z; ... chopped sum 
355 
y0:=y00x00100000; ... y := y/2 is correctly rounded. 
356 
I := i; ... restore inexact flag 
357 
R := r; ... restore rounded mode 
358 
return sqrt(x):=y. 
359 

360 
(4) Special cases 
361 

362 
Square root of +inf, +0, or NaN is itself; 
363 
Square root of a negative number is NaN with invalid signal. 
364 

365 

366 
B. sqrt(x) by Reciproot Iteration 
367 

368 
(1) Initial approximation 
369 

370 
Let x0 and x1 be the leading and the trailing 32bit words of 
371 
a floating point number x (in IEEE double format) respectively 
372 
(see section A). By performing shifs and subtracts on x0 and y0, 
373 
we obtain a 7.8bit approximation of 1/sqrt(x) as follows. 
374 

375 
k := 0x5fe80000  (x0>>1); 
376 
y0:= k  T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits 
377 

378 
Here k is a 32bit integer and T2[] is an integer array 
379 
containing correction terms. Now magically the floating 
380 
value of y (y's leading 32bit word is y0, the value of 
381 
its trailing word y1 is set to zero) approximates 1/sqrt(x) 
382 
to almost 7.8bit. 
383 

384 
Value of T2: 
385 
static int T2[64]= { 
386 
0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866, 
387 
0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f, 
388 
0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d, 
389 
0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0, 
390 
0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989, 
391 
0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd, 
392 
0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e, 
393 
0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,}; 
394 

395 
(2) Iterative refinement 
396 

397 
Apply Reciproot iteration three times to y and multiply the 
398 
result by x to get an approximation z that matches sqrt(x) 
399 
to about 1 ulp. To be exact, we will have 
400 
1ulp < sqrt(x)z<1.0625ulp. 
401 

402 
... set rounding mode to Roundtonearest 
403 
y := y*(1.50.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x) 
404 
y := y*((1.52^30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x) 
405 
... special arrangement for better accuracy 
406 
z := x*y ... 29 bits to sqrt(x), with z*y<1 
407 
z := z + 0.5*z*(1z*y) ... about 1 ulp to sqrt(x) 
408 

409 
Remark 2. The constant 1.52^30 is chosen to bias the error so that 
410 
(a) the term z*y in the final iteration is always less than 1; 
411 
(b) the error in the final result is biased upward so that 
412 
1 ulp < sqrt(x)  z < 1.0625 ulp 
413 
instead of sqrt(x)z<1.03125ulp. 
414 

415 
(3) Final adjustment 
416 

417 
By twiddling y's last bit it is possible to force y to be 
418 
correctly rounded according to the prevailing rounding mode 
419 
as follows. Let r and i be copies of the rounding mode and 
420 
inexact flag before entering the square root program. Also we 
421 
use the expression y+ulp for the next representable floating 
422 
numbers (up and down) of y. Note that y+ulp = either fixed 
423 
point y+1, or multiply y by nextafter(1,+inf) in chopped 
424 
mode. 
425 

426 
R := RZ; ... set rounding mode to roundtowardzero 
427 
switch(r) { 
428 
case RN: ... roundtonearest 
429 
if(x<= z*(zulp)...chopped) z = z  ulp; else 
430 
if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp; 
431 
break; 
432 
case RZ:case RM: ... roundtozero or roundtoinf 
433 
R:=RP; ... reset rounding mod to roundto+inf 
434 
if(x<z*z ... rounded up) z = z  ulp; else 
435 
if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp; 
436 
break; 
437 
case RP: ... roundto+inf 
438 
if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else 
439 
if(x>z*z ...chopped) z = z+ulp; 
440 
break; 
441 
} 
442 

443 
Remark 3. The above comparisons can be done in fixed point. For 
444 
example, to compare x and w=z*z chopped, it suffices to compare 
445 
x1 and w1 (the trailing parts of x and w), regarding them as 
446 
two's complement integers. 
447 

448 
...Is z an exact square root? 
449 
To determine whether z is an exact square root of x, let z1 be the 
450 
trailing part of z, and also let x0 and x1 be the leading and 
451 
trailing parts of x. 
452 

453 
If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0 
454 
I := 1; ... Raise Inexact flag: z is not exact 
455 
else { 
456 
j := 1  [(x0>>20)&1] ... j = logb(x) mod 2 
457 
k := z1 >> 26; ... get z's 25th and 26th 
458 
fraction bits 
459 
I := i or (k&j) or ((k&(j+j+1))!=(x1&3)); 
460 
} 
461 
R:= r ... restore rounded mode 
462 
return sqrt(x):=z. 
463 

464 
If multiplication is cheaper then the foregoing red tape, the 
465 
Inexact flag can be evaluated by 
466 

467 
I := i; 
468 
I := (z*z!=x) or I. 
469 

470 
Note that z*z can overwrite I; this value must be sensed if it is 
471 
True. 
472 

473 
Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be 
474 
zero. 
475 

476 
 
477 
z1:  f2  
478 
 
479 
bit 31 bit 0 
480 

481 
Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd 
482 
or even of logb(x) have the following relations: 
483 

484 
 
485 
bit 27,26 of z1 bit 1,0 of x1 logb(x) 
486 
 
487 
00 00 odd and even 
488 
01 01 even 
489 
10 10 odd 
490 
10 00 even 
491 
11 01 even 
492 
 
493 

494 
(4) Special cases (see (4) of Section A). 
495 

496 
*/ 
497 
