src/libm/e_sqrt.c
author Ryan C. Gordon <icculus@icculus.org>
Mon, 21 May 2018 12:00:21 -0400
changeset 11994 8e094f91ddab
parent 11683 48bcba563d9c
permissions -rw-r--r--
thread: fixed compiler warnings on non-Linux systems that use pthread.

(static function rtkit_setpriority was unused, moved it in with rest of
__LINUX__ section.)
     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[32]= {
   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[64]= {
   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  */