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