sdl-ios-xcode: src/libm/e

annotate src/libm/e_sqrt.c @ 3294:470d0a416aa7

Fixed bug #714 fuzzyTew@gmail.com 2009-03-14 15:18:45 PDT patch to change HAVE_ICONV to HAVE_ICONV_H There are two separate iconv checks in configure -- one for the header file and one for the library. include/SDL_stdinc.h uses the library define to see whether or not it should reference the types defined in the header, which naturally breaks if the library exists and the header does not.

author	Sam Lantinga <slouken@libsdl.org>
date	Mon, 21 Sep 2009 11:04:01 +0000
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 = 2q , 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.5xyy) ... 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.5xyy)... 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 := xy ... 29 bits to sqrt(x), with zy<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.5z(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+2ulp; 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 */

Mercurial > sdl-ios-xcode

annotate src/libm/e_sqrt.c @ 3294:470d0a416aa7