Mercurial > sdl-ios-xcode

annotate src/libm/e_sqrt.c @ 4573:6399178be313

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Completed work on X11_CreateTexture. Added lots of safety features. These include support for drawing a texture using the core protocol while other textures are drawn using Xrender if Xrender does not support the color format of the said texture or any other fault with Xrender.
author Sunny Sachanandani <sunnysachanandani@gmail.com>
date Fri, 28 May 2010 20:40:09 +0530
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 */