libm/src/e_sqrt.c

13 /* sqrt(x)
14  * Return correctly rounded sqrt.
16  *	     |  Use the hardware sqrt if you have one |
23  *		sqrt(x) = 2^k * sqrt(y)
25  *	Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
74  *	sqrt(+-0) = +-0 	... exact
75  *	sqrt(inf) = inf
76  *	sqrt(-ve) = NaN		... with invalid signal
77  *	sqrt(NaN) = NaN		... with invalid signal for signaling NaN
91 sqrt(double x)
102 	    return x*x+x;		/* sqrt(NaN)=NaN, sqrt(+inf)=+inf
103 					   sqrt(-inf)=sNaN */
107 	    if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
109 		return (x-x)/(x-x);		/* sqrt(-ve) = sNaN */
131     /* generate sqrt(x) bit by bit */
134 	q = q1 = s0 = s1 = 0;	/* [q,q1] = sqrt(x) */
186 DEF_STD(sqrt);
187 LDBL_MAYBE_CLONE(sqrt);
195 	Two algorithms are given here to implement sqrt(x)
197 	Both supply sqrt(x) correctly rounded. The first algorithm (in
208 A.  sqrt(x) by Newton Iteration
227 	as integers), we obtain an 8-bit approximation of sqrt(x) as
231 		y0 := k - T1[31&(k>>15)].	... y ~ sqrt(x) to 8 bits
235 	approximates sqrt(x) to almost 8-bit.
247 	sqrt(x) to within 1 ulp (Unit in the Last Place):
257 		y := (y+x/y)		... almost 17 sig. bits to 2*sqrt(x)
258 		y := y - 0x00100006	... almost 18 sig. bits to sqrt(x)
292 		        return sqrt(x):=y.
297 		i := TRUE;		... sqrt(x) is inexact
306 	        return sqrt(x):=y.
314 B.  sqrt(x) by Reciproot Iteration
321 	we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
324 	    y0:= k - T2[63&(k>>14)].	... y ~ 1/sqrt(x) to 7.8 bits
329 	its trailing word y1 is set to zero) approximates 1/sqrt(x)
346 	result by x to get an approximation z that matches sqrt(x)
348 		-1ulp < sqrt(x)-z<1.0625ulp.
351 	   y := y*(1.5-0.5*x*y*y)	... almost 15 sig. bits to 1/sqrt(x)
352 	   y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
354 	   z := x*y			... 29 bits to sqrt(x), with z*y<1
355 	   z := z + 0.5*z*(1-z*y)	... about 1 ulp to sqrt(x)
360 		-1 ulp < sqrt(x) - z < 1.0625 ulp
361 	    instead of |sqrt(x)-z|<1.03125ulp.
410 	return sqrt(x):=z.