Mercurial > sdl-ios-xcode

comparison src/video/e_sqrt.h @ 1330:450721ad5436

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