sdl-ios-xcode: src/libm/e

annotate src/libm/e_sqrt.c @ 2842:97ba0be8b565

Date: Sat, 06 Dec 2008 15:27:00 +0100 From: Couriersud Subject: SDL: Relative mouse movements The patch below will reenable processing of relative mouse movements. The DirectFB drivers generates those in "grabbed" mode. These ensure, that even in fullscreen mode relative movements are reported. SDLMAME depends on this for games with trackballs. Looking at the code I ask myself whether relative movements should be handled in the drivers (x11, directfb). Both x11 and directfb are able to report relative movements. This would leave it to the driver to use the most appropriate method for relative movements when at the border of a fullscreen window or being "grabbed".

author	Sam Lantinga <slouken@libsdl.org>
date	Sat, 06 Dec 2008 17:50:50 +0000
parents	a98604b691c8
children	dc1eb82ffdaa

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)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	14 static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	15 #endif
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	16
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	17 /* __ieee754_sqrt(x)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	18 * Return correctly rounded sqrt.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	19 * ------------------------------------------
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	20 * \| 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	21 * ------------------------------------------
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	22 * Method:
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	23 * 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	24 * 1. Normalization
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	25 * 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	26 * 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	27 * 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	28 * 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	29 * 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	30 * i 0
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	31 * i+1 2
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	32 * 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	33 * i i i i
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	34 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	35 * 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	36 * i+1 i
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	37 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	38 * -(i+1) 2
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	39 * (q + 2 ) <= y. (2)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	40 * i
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	41 * -(i+1)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	42 * 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	43 * 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	44 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	45 * 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	46 * 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	47 * -(i+1)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	48 * s + 2 <= y (3)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	49 * i i
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	50 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	51 * 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	52 * i i
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	53 * the following recurrence formula:
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	54 * if (3) is false
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	55 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	56 * 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	57 * 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	58 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	59 * otherwise,
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	60 * -i -(i+1)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	61 * 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	62 * 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	63 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	64 * 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	65 * 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	66 * 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	67 * in (3).
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	68 * 3. Final rounding
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	69 * 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	70 * 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	71 * 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	72 * (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	73 * 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	74 * 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	75 * 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	76 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	77 * Special cases:
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	78 * sqrt(+-0) = +-0 ... exact
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	79 * sqrt(inf) = inf
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	80 * 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	81 * 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	82 *
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	83 * 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	84 *---------------
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 #include "math.h"
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	88 #include "math_private.h"
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	89
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	90 #ifdef __STDC__
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	91 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	92 #else
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	93 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	94 #endif
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	95
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	96 #ifdef __STDC__
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	97 double attribute_hidden
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	98 __ieee754_sqrt(double x)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	99 #else
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	100 double attribute_hidden
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	101 __ieee754_sqrt(x)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	102 double x;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	103 #endif
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	104 {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	105 double z;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	106 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	107 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	108 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	109
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	110 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	111
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	112 /* 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	113 if ((ix0 & 0x7ff00000) == 0x7ff00000) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	114 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	115 sqrt(-inf)=sNaN */
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	116 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	117 /* take care of zero */
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	118 if (ix0 <= 0) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	119 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	120 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	121 else if (ix0 < 0)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	122 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	123 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	124 /* normalize x */
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	125 m = (ix0 >> 20);
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	126 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	127 while (ix0 == 0) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	128 m -= 21;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	129 ix0 \|= (ix1 >> 11);
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	130 ix1 <<= 21;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	131 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	132 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	133 ix0 <<= 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	134 m -= i - 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	135 ix0 \|= (ix1 >> (32 - i));
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	136 ix1 <<= i;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	137 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	138 m -= 1023; /* unbias exponent */
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	139 ix0 = (ix0 & 0x000fffff) \| 0x00100000;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	140 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	141 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	142 ix1 += ix1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	143 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	144 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	145
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	146 /* 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	147 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	148 ix1 += ix1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	149 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	150 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	151
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	152 while (r != 0) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	153 t = s0 + r;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	154 if (t <= ix0) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	155 s0 = t + r;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	156 ix0 -= t;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	157 q += r;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	158 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	159 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	160 ix1 += ix1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	161 r >>= 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	162 }
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 r = sign;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	165 while (r != 0) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	166 t1 = s1 + r;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	167 t = s0;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	168 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	169 s1 = t1 + r;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	170 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	171 s0 += 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	172 ix0 -= t;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	173 if (ix1 < t1)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	174 ix0 -= 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	175 ix1 -= t1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	176 q1 += r;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	177 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	178 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	179 ix1 += ix1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	180 r >>= 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	181 }
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 /* 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	184 if ((ix0 \| ix1) != 0) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	185 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	186 if (z >= one) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	187 z = one + tiny;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	188 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	189 q1 = 0;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	190 q += 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	191 } else if (z > one) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	192 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	193 q += 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	194 q1 += 2;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	195 } else
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	196 q1 += (q1 & 1);
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	197 }
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 ix0 = (q >> 1) + 0x3fe00000;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	200 ix1 = q1 >> 1;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	201 if ((q & 1) == 1)
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	202 ix1 \|= sign;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	203 ix0 += (m << 20);
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	204 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	205 return z;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	206 }
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 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	210 -------------
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	211 (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	212 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	213
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	214 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	215 (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	216 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	217 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	218 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	219 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	220 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	221 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	222 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	223 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	224 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	225 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	226
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	227 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	228
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	229 (1) Initial approximation
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	230
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	231 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	232 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	233
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	234 1 11 52 ...widths
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	235 ------------------------------------------------------
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	236 x: \|s\| e \| f \|
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	237 ------------------------------------------------------
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	238 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	239
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 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	243 ------------------------ ------------------------
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 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	246 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	247 follows.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	248
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	249 k := (x0>>1) + 0x1ff80000;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	250 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	251 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	252 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	253 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	254 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	255
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	256 Value of T1:
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	257 static int T1[32]= {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	258 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	259 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	260 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	261 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	262
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	263 (2) Iterative refinement
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	264
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	265 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	266 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	267
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	268 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	269 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	270 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	271
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 Remark 1.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	274 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	275
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	276 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	277 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	278
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	279 2
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	280 (x-y )*y
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	281 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	282 2
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	283 3y + x
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	284
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 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	287 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	288 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	289 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	290 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	291 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	292
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	293 (3) Final adjustment
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	294
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	295 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	296 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	297 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	298 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	299 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	300 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	301 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	302 mode.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	303
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	304 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	305 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	306 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	307 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	308 if(z=y) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	309 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	310 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	311 return sqrt(x):=y.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	312 } else {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	313 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	314 }
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 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	317 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	318 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	319 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	320 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	321 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	322 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	323 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	324 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	325 return sqrt(x):=y.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	326
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	327 (4) Special cases
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	328
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	329 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	330 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	331
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 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	334
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	335 (1) Initial approximation
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	336
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	337 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	338 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	339 (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	340 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	341
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	342 k := 0x5fe80000 - (x0>>1);
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	343 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	344
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	345 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	346 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	347 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	348 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	349 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	350
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	351 Value of T2:
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	352 static int T2[64]= {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	353 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	354 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	355 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	356 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	357 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	358 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	359 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	360 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	361
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	362 (2) Iterative refinement
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	363
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	364 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	365 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	366 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	367 -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	368
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	369 ... 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	370 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	371 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	372 ... 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	373 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	374 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	375
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	376 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	377 (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	378 (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	379 -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	380 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	381
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	382 (3) Final adjustment
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	383
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	384 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	385 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	386 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	387 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	388 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	389 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	390 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	391 mode.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	392
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	393 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	394 switch(r) {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	395 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	396 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	397 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	398 break;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	399 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	400 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	401 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	402 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	403 break;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	404 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	405 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	406 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	407 break;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	408 }
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 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	411 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	412 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	413 two's complement integers.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	414
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	415 ...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	416 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	417 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	418 trailing parts of x.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	419
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	420 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	421 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	422 else {
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	423 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	424 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	425 fraction bits
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	426 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	427 }
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	428 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	429 return sqrt(x):=z.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	430
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	431 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	432 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	433
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	434 I := i;
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	435 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	436
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	437 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	438 True.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	439
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	440 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	441 zero.
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	442
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 z1: \| f2 \|
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	445 --------------------
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	446 bit 31 bit 0
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	447
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	448 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	449 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	450
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 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	453 -------------------------------------------------
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	454 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	455 01 01 even
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	456 10 10 odd
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	457 10 00 even
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	458 11 01 even
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	459 -------------------------------------------------
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 (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	462
a98604b691c8 Expanded the libm support and put it into a separate directory. Sam Lantinga <slouken@libsdl.org> parents: diff changeset	463 */

Mercurial > sdl-ios-xcode

annotate src/libm/e_sqrt.c @ 2842:97ba0be8b565