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