1 /** @file
2 Compute the logrithm of x.
3
4 Copyright (c) 2010 - 2011, Intel Corporation. All rights reserved.<BR>
5 This program and the accompanying materials are licensed and made available under
6 the terms and conditions of the BSD License that accompanies this distribution.
7 The full text of the license may be found at
8 http://opensource.org/licenses/bsd-license.
9
10 THE PROGRAM IS DISTRIBUTED UNDER THE BSD LICENSE ON AN "AS IS" BASIS,
11 WITHOUT WARRANTIES OR REPRESENTATIONS OF ANY KIND, EITHER EXPRESS OR IMPLIED.
12
13 * ====================================================
14 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
15 *
16 * Developed at SunPro, a Sun Microsystems, Inc. business.
17 * Permission to use, copy, modify, and distribute this
18 * software is freely granted, provided that this notice
19 * is preserved.
20 * ====================================================
21
22 e_sqrt.c 5.1 93/09/24
23 NetBSD: e_sqrt.c,v 1.12 2002/05/26 22:01:52 wiz Exp
24 **/
25 #include <LibConfig.h>
26 #include <sys/EfiCdefs.h>
27
28 #include <errno.h>
29 #include "math.h"
30 #include "math_private.h"
31
32 #if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */
33 // potential divide by 0 -- near line 129, (x-x)/(x-x) is on purpose
34 #pragma warning ( disable : 4723 )
35 #endif
36
37 /* __ieee754_sqrt(x)
38 * Return correctly rounded sqrt.
39 * ------------------------------------------
40 * | Use the hardware sqrt if you have one |
41 * ------------------------------------------
42 * Method:
43 * Bit by bit method using integer arithmetic. (Slow, but portable)
44 * 1. Normalization
45 * Scale x to y in [1,4) with even powers of 2:
46 * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
47 * sqrt(x) = 2^k * sqrt(y)
48 * 2. Bit by bit computation
49 * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
50 * i 0
51 * i+1 2
52 * s = 2*q , and y = 2 * ( y - q ). (1)
53 * i i i i
54 *
55 * To compute q from q , one checks whether
56 * i+1 i
57 *
58 * -(i+1) 2
59 * (q + 2 ) <= y. (2)
60 * i
61 * -(i+1)
62 * If (2) is false, then q = q ; otherwise q = q + 2 .
63 * i+1 i i+1 i
64 *
65 * With some algebric manipulation, it is not difficult to see
66 * that (2) is equivalent to
67 * -(i+1)
68 * s + 2 <= y (3)
69 * i i
70 *
71 * The advantage of (3) is that s and y can be computed by
72 * i i
73 * the following recurrence formula:
74 * if (3) is false
75 *
76 * s = s , y = y ; (4)
77 * i+1 i i+1 i
78 *
79 * otherwise,
80 * -i -(i+1)
81 * s = s + 2 , y = y - s - 2 (5)
82 * i+1 i i+1 i i
83 *
84 * One may easily use induction to prove (4) and (5).
85 * Note. Since the left hand side of (3) contain only i+2 bits,
86 * it does not necessary to do a full (53-bit) comparison
87 * in (3).
88 * 3. Final rounding
89 * After generating the 53 bits result, we compute one more bit.
90 * Together with the remainder, we can decide whether the
91 * result is exact, bigger than 1/2ulp, or less than 1/2ulp
92 * (it will never equal to 1/2ulp).
93 * The rounding mode can be detected by checking whether
94 * huge + tiny is equal to huge, and whether huge - tiny is
95 * equal to huge for some floating point number "huge" and "tiny".
96 *
97 * Special cases:
98 * sqrt(+-0) = +-0 ... exact
99 * sqrt(inf) = inf
100 * sqrt(-ve) = NaN ... with invalid signal
101 * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
102 *
103 * Other methods : see the appended file at the end of the program below.
104 *---------------
105 */
106
107 static const double one = 1.0, tiny=1.0e-300;
108
109 double
__ieee754_sqrt(double x)110 __ieee754_sqrt(double x)
111 {
112 double z;
113 int32_t sign = (int)0x80000000;
114 int32_t ix0,s0,q,m,t,i;
115 u_int32_t r,t1,s1,ix1,q1;
116
117 EXTRACT_WORDS(ix0,ix1,x);
118
119 /* take care of Inf and NaN */
120 if((ix0&0x7ff00000)==0x7ff00000) {
121 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
122 sqrt(-inf)=sNaN */
123 }
124 /* take care of zero */
125 if(ix0<=0) {
126 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
127 else if(ix0<0) {
128 errno = EDOM;
129 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
130 }
131 }
132 /* normalize x */
133 m = (ix0>>20);
134 if(m==0) { /* subnormal x */
135 while(ix0==0) {
136 m -= 21;
137 ix0 |= (ix1>>11); ix1 <<= 21;
138 }
139 for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
140 m -= i-1;
141 ix0 |= (ix1>>(32-i));
142 ix1 <<= i;
143 }
144 m -= 1023; /* unbias exponent */
145 ix0 = (ix0&0x000fffff)|0x00100000;
146 if(m&1){ /* odd m, double x to make it even */
147 ix0 += ix0 + ((ix1&sign)>>31);
148 ix1 += ix1;
149 }
150 m >>= 1; /* m = [m/2] */
151
152 /* generate sqrt(x) bit by bit */
153 ix0 += ix0 + ((ix1&sign)>>31);
154 ix1 += ix1;
155 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
156 r = 0x00200000; /* r = moving bit from right to left */
157
158 while(r!=0) {
159 t = s0+r;
160 if(t<=ix0) {
161 s0 = t+r;
162 ix0 -= t;
163 q += r;
164 }
165 ix0 += ix0 + ((ix1&sign)>>31);
166 ix1 += ix1;
167 r>>=1;
168 }
169
170 r = sign;
171 while(r!=0) {
172 t1 = s1+r;
173 t = s0;
174 if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
175 s1 = t1+r;
176 if(((t1&sign)==(u_int32_t)sign)&&(s1&sign)==0) s0 += 1;
177 ix0 -= t;
178 if (ix1 < t1) ix0 -= 1;
179 ix1 -= t1;
180 q1 += r;
181 }
182 ix0 += ix0 + ((ix1&sign)>>31);
183 ix1 += ix1;
184 r>>=1;
185 }
186
187 /* use floating add to find out rounding direction */
188 if((ix0|ix1)!=0) {
189 z = one-tiny; /* trigger inexact flag */
190 if (z>=one) {
191 z = one+tiny;
192 if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
193 else if (z>one) {
194 if (q1==(u_int32_t)0xfffffffe) q+=1;
195 q1+=2;
196 } else
197 q1 += (q1&1);
198 }
199 }
200 ix0 = (q>>1)+0x3fe00000;
201 ix1 = q1>>1;
202 if ((q&1)==1) ix1 |= sign;
203 ix0 += (m <<20);
204 INSERT_WORDS(z,ix0,ix1);
205 return z;
206 }
207
208 /*
209 Other methods (use floating-point arithmetic)
210 -------------
211 (This is a copy of a drafted paper by Prof W. Kahan
212 and K.C. Ng, written in May, 1986)
213
214 Two algorithms are given here to implement sqrt(x)
215 (IEEE double precision arithmetic) in software.
216 Both supply sqrt(x) correctly rounded. The first algorithm (in
217 Section A) uses newton iterations and involves four divisions.
218 The second one uses reciproot iterations to avoid division, but
219 requires more multiplications. Both algorithms need the ability
220 to chop results of arithmetic operations instead of round them,
221 and the INEXACT flag to indicate when an arithmetic operation
222 is executed exactly with no roundoff error, all part of the
223 standard (IEEE 754-1985). The ability to perform shift, add,
224 subtract and logical AND operations upon 32-bit words is needed
225 too, though not part of the standard.
226
227 A. sqrt(x) by Newton Iteration
228
229 (1) Initial approximation
230
231 Let x0 and x1 be the leading and the trailing 32-bit words of
232 a floating point number x (in IEEE double format) respectively
233
234 1 11 52 ...widths
235 ------------------------------------------------------
236 x: |s| e | f |
237 ------------------------------------------------------
238 msb lsb msb lsb ...order
239
240
241 ------------------------ ------------------------
242 x0: |s| e | f1 | x1: | f2 |
243 ------------------------ ------------------------
244
245 By performing shifts and subtracts on x0 and x1 (both regarded
246 as integers), we obtain an 8-bit approximation of sqrt(x) as
247 follows.
248
249 k := (x0>>1) + 0x1ff80000;
250 y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
251 Here k is a 32-bit integer and T1[] is an integer array containing
252 correction terms. Now magically the floating value of y (y's
253 leading 32-bit word is y0, the value of its trailing word is 0)
254 approximates sqrt(x) to almost 8-bit.
255
256 Value of T1:
257 static int T1[32]= {
258 0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
259 29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
260 83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
261 16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
262
263 (2) Iterative refinement
264
265 Apply Heron's rule three times to y, we have y approximates
266 sqrt(x) to within 1 ulp (Unit in the Last Place):
267
268 y := (y+x/y)/2 ... almost 17 sig. bits
269 y := (y+x/y)/2 ... almost 35 sig. bits
270 y := y-(y-x/y)/2 ... within 1 ulp
271
272
273 Remark 1.
274 Another way to improve y to within 1 ulp is:
275
276 y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
277 y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
278
279 2
280 (x-y )*y
281 y := y + 2* ---------- ...within 1 ulp
282 2
283 3y + x
284
285
286 This formula has one division fewer than the one above; however,
287 it requires more multiplications and additions. Also x must be
288 scaled in advance to avoid spurious overflow in evaluating the
289 expression 3y*y+x. Hence it is not recommended uless division
290 is slow. If division is very slow, then one should use the
291 reciproot algorithm given in section B.
292
293 (3) Final adjustment
294
295 By twiddling y's last bit it is possible to force y to be
296 correctly rounded according to the prevailing rounding mode
297 as follows. Let r and i be copies of the rounding mode and
298 inexact flag before entering the square root program. Also we
299 use the expression y+-ulp for the next representable floating
300 numbers (up and down) of y. Note that y+-ulp = either fixed
301 point y+-1, or multiply y by nextafter(1,+-inf) in chopped
302 mode.
303
304 I := FALSE; ... reset INEXACT flag I
305 R := RZ; ... set rounding mode to round-toward-zero
306 z := x/y; ... chopped quotient, possibly inexact
307 If(not I) then { ... if the quotient is exact
308 if(z=y) {
309 I := i; ... restore inexact flag
310 R := r; ... restore rounded mode
311 return sqrt(x):=y.
312 } else {
313 z := z - ulp; ... special rounding
314 }
315 }
316 i := TRUE; ... sqrt(x) is inexact
317 If (r=RN) then z=z+ulp ... rounded-to-nearest
318 If (r=RP) then { ... round-toward-+inf
319 y = y+ulp; z=z+ulp;
320 }
321 y := y+z; ... chopped sum
322 y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
323 I := i; ... restore inexact flag
324 R := r; ... restore rounded mode
325 return sqrt(x):=y.
326
327 (4) Special cases
328
329 Square root of +inf, +-0, or NaN is itself;
330 Square root of a negative number is NaN with invalid signal.
331
332
333 B. sqrt(x) by Reciproot Iteration
334
335 (1) Initial approximation
336
337 Let x0 and x1 be the leading and the trailing 32-bit words of
338 a floating point number x (in IEEE double format) respectively
339 (see section A). By performing shifs and subtracts on x0 and y0,
340 we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
341
342 k := 0x5fe80000 - (x0>>1);
343 y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
344
345 Here k is a 32-bit integer and T2[] is an integer array
346 containing correction terms. Now magically the floating
347 value of y (y's leading 32-bit word is y0, the value of
348 its trailing word y1 is set to zero) approximates 1/sqrt(x)
349 to almost 7.8-bit.
350
351 Value of T2:
352 static int T2[64]= {
353 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
354 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
355 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
356 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
357 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
358 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
359 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
360 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
361
362 (2) Iterative refinement
363
364 Apply Reciproot iteration three times to y and multiply the
365 result by x to get an approximation z that matches sqrt(x)
366 to about 1 ulp. To be exact, we will have
367 -1ulp < sqrt(x)-z<1.0625ulp.
368
369 ... set rounding mode to Round-to-nearest
370 y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
371 y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
372 ... special arrangement for better accuracy
373 z := x*y ... 29 bits to sqrt(x), with z*y<1
374 z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
375
376 Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
377 (a) the term z*y in the final iteration is always less than 1;
378 (b) the error in the final result is biased upward so that
379 -1 ulp < sqrt(x) - z < 1.0625 ulp
380 instead of |sqrt(x)-z|<1.03125ulp.
381
382 (3) Final adjustment
383
384 By twiddling y's last bit it is possible to force y to be
385 correctly rounded according to the prevailing rounding mode
386 as follows. Let r and i be copies of the rounding mode and
387 inexact flag before entering the square root program. Also we
388 use the expression y+-ulp for the next representable floating
389 numbers (up and down) of y. Note that y+-ulp = either fixed
390 point y+-1, or multiply y by nextafter(1,+-inf) in chopped
391 mode.
392
393 R := RZ; ... set rounding mode to round-toward-zero
394 switch(r) {
395 case RN: ... round-to-nearest
396 if(x<= z*(z-ulp)...chopped) z = z - ulp; else
397 if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
398 break;
399 case RZ:case RM: ... round-to-zero or round-to--inf
400 R:=RP; ... reset rounding mod to round-to-+inf
401 if(x<z*z ... rounded up) z = z - ulp; else
402 if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
403 break;
404 case RP: ... round-to-+inf
405 if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
406 if(x>z*z ...chopped) z = z+ulp;
407 break;
408 }
409
410 Remark 3. The above comparisons can be done in fixed point. For
411 example, to compare x and w=z*z chopped, it suffices to compare
412 x1 and w1 (the trailing parts of x and w), regarding them as
413 two's complement integers.
414
415 ...Is z an exact square root?
416 To determine whether z is an exact square root of x, let z1 be the
417 trailing part of z, and also let x0 and x1 be the leading and
418 trailing parts of x.
419
420 If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
421 I := 1; ... Raise Inexact flag: z is not exact
422 else {
423 j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
424 k := z1 >> 26; ... get z's 25-th and 26-th
425 fraction bits
426 I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
427 }
428 R:= r ... restore rounded mode
429 return sqrt(x):=z.
430
431 If multiplication is cheaper than the foregoing red tape, the
432 Inexact flag can be evaluated by
433
434 I := i;
435 I := (z*z!=x) or I.
436
437 Note that z*z can overwrite I; this value must be sensed if it is
438 True.
439
440 Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
441 zero.
442
443 --------------------
444 z1: | f2 |
445 --------------------
446 bit 31 bit 0
447
448 Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
449 or even of logb(x) have the following relations:
450
451 -------------------------------------------------
452 bit 27,26 of z1 bit 1,0 of x1 logb(x)
453 -------------------------------------------------
454 00 00 odd and even
455 01 01 even
456 10 10 odd
457 10 00 even
458 11 01 even
459 -------------------------------------------------
460
461 (4) Special cases (see (4) of Section A).
462
463 */
464
465