annotate src/libm/e_sqrt.c @ 5221:16d5e8e90e01

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