src/libm/e_sqrt.c
changeset 11683 48bcba563d9c
parent 6044 35448a5ea044
     1.1 --- a/src/libm/e_sqrt.c	Sat Nov 04 15:34:14 2017 -0700
     1.2 +++ b/src/libm/e_sqrt.c	Sat Nov 04 15:53:19 2017 -0700
     1.3 @@ -1,4 +1,3 @@
     1.4 -/* @(#)e_sqrt.c 5.1 93/09/24 */
     1.5  /*
     1.6   * ====================================================
     1.7   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     1.8 @@ -10,11 +9,6 @@
     1.9   * ====================================================
    1.10   */
    1.11  
    1.12 -#if defined(LIBM_SCCS) && !defined(lint)
    1.13 -static const char rcsid[] =
    1.14 -    "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
    1.15 -#endif
    1.16 -
    1.17  /* __ieee754_sqrt(x)
    1.18   * Return correctly rounded sqrt.
    1.19   *           ------------------------------------------
    1.20 @@ -88,125 +82,124 @@
    1.21  #include "math_libm.h"
    1.22  #include "math_private.h"
    1.23  
    1.24 -#ifdef __STDC__
    1.25  static const double one = 1.0, tiny = 1.0e-300;
    1.26 -#else
    1.27 -static double one = 1.0, tiny = 1.0e-300;
    1.28 -#endif
    1.29  
    1.30 -#ifdef __STDC__
    1.31 -double attribute_hidden
    1.32 -__ieee754_sqrt(double x)
    1.33 -#else
    1.34 -double attribute_hidden
    1.35 -__ieee754_sqrt(x)
    1.36 -     double x;
    1.37 -#endif
    1.38 +double attribute_hidden __ieee754_sqrt(double x)
    1.39  {
    1.40 -    double z;
    1.41 -    int32_t sign = (int) 0x80000000;
    1.42 -    int32_t ix0, s0, q, m, t, i;
    1.43 -    u_int32_t r, t1, s1, ix1, q1;
    1.44 +	double z;
    1.45 +	int32_t sign = (int)0x80000000;
    1.46 +	int32_t ix0,s0,q,m,t,i;
    1.47 +	u_int32_t r,t1,s1,ix1,q1;
    1.48  
    1.49 -    EXTRACT_WORDS(ix0, ix1, x);
    1.50 +	EXTRACT_WORDS(ix0,ix1,x);
    1.51  
    1.52      /* take care of Inf and NaN */
    1.53 -    if ((ix0 & 0x7ff00000) == 0x7ff00000) {
    1.54 -        return x * x + x;       /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
    1.55 -                                   sqrt(-inf)=sNaN */
    1.56 -    }
    1.57 +	if((ix0&0x7ff00000)==0x7ff00000) {
    1.58 +	    return x*x+x;		/* sqrt(NaN)=NaN, sqrt(+inf)=+inf
    1.59 +					   sqrt(-inf)=sNaN */
    1.60 +	}
    1.61      /* take care of zero */
    1.62 -    if (ix0 <= 0) {
    1.63 -        if (((ix0 & (~sign)) | ix1) == 0)
    1.64 -            return x;           /* sqrt(+-0) = +-0 */
    1.65 -        else if (ix0 < 0)
    1.66 -            return (x - x) / (x - x);   /* sqrt(-ve) = sNaN */
    1.67 -    }
    1.68 +	if(ix0<=0) {
    1.69 +	    if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
    1.70 +	    else if(ix0<0)
    1.71 +		return (x-x)/(x-x);		/* sqrt(-ve) = sNaN */
    1.72 +	}
    1.73      /* normalize x */
    1.74 -    m = (ix0 >> 20);
    1.75 -    if (m == 0) {               /* subnormal x */
    1.76 -        while (ix0 == 0) {
    1.77 -            m -= 21;
    1.78 -            ix0 |= (ix1 >> 11);
    1.79 -            ix1 <<= 21;
    1.80 -        }
    1.81 -        for (i = 0; (ix0 & 0x00100000) == 0; i++)
    1.82 -            ix0 <<= 1;
    1.83 -        m -= i - 1;
    1.84 -        ix0 |= (ix1 >> (32 - i));
    1.85 -        ix1 <<= i;
    1.86 -    }
    1.87 -    m -= 1023;                  /* unbias exponent */
    1.88 -    ix0 = (ix0 & 0x000fffff) | 0x00100000;
    1.89 -    if (m & 1) {                /* odd m, double x to make it even */
    1.90 -        ix0 += ix0 + ((ix1 & sign) >> 31);
    1.91 -        ix1 += ix1;
    1.92 -    }
    1.93 -    m >>= 1;                    /* m = [m/2] */
    1.94 +	m = (ix0>>20);
    1.95 +	if(m==0) {				/* subnormal x */
    1.96 +	    while(ix0==0) {
    1.97 +		m -= 21;
    1.98 +		ix0 |= (ix1>>11); ix1 <<= 21;
    1.99 +	    }
   1.100 +	    for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
   1.101 +	    m -= i-1;
   1.102 +	    ix0 |= (ix1>>(32-i));
   1.103 +	    ix1 <<= i;
   1.104 +	}
   1.105 +	m -= 1023;	/* unbias exponent */
   1.106 +	ix0 = (ix0&0x000fffff)|0x00100000;
   1.107 +	if(m&1){	/* odd m, double x to make it even */
   1.108 +	    ix0 += ix0 + ((ix1&sign)>>31);
   1.109 +	    ix1 += ix1;
   1.110 +	}
   1.111 +	m >>= 1;	/* m = [m/2] */
   1.112  
   1.113      /* generate sqrt(x) bit by bit */
   1.114 -    ix0 += ix0 + ((ix1 & sign) >> 31);
   1.115 -    ix1 += ix1;
   1.116 -    q = q1 = s0 = s1 = 0;       /* [q,q1] = sqrt(x) */
   1.117 -    r = 0x00200000;             /* r = moving bit from right to left */
   1.118 +	ix0 += ix0 + ((ix1&sign)>>31);
   1.119 +	ix1 += ix1;
   1.120 +	q = q1 = s0 = s1 = 0;	/* [q,q1] = sqrt(x) */
   1.121 +	r = 0x00200000;		/* r = moving bit from right to left */
   1.122  
   1.123 -    while (r != 0) {
   1.124 -        t = s0 + r;
   1.125 -        if (t <= ix0) {
   1.126 -            s0 = t + r;
   1.127 -            ix0 -= t;
   1.128 -            q += r;
   1.129 -        }
   1.130 -        ix0 += ix0 + ((ix1 & sign) >> 31);
   1.131 -        ix1 += ix1;
   1.132 -        r >>= 1;
   1.133 -    }
   1.134 +	while(r!=0) {
   1.135 +	    t = s0+r;
   1.136 +	    if(t<=ix0) {
   1.137 +		s0   = t+r;
   1.138 +		ix0 -= t;
   1.139 +		q   += r;
   1.140 +	    }
   1.141 +	    ix0 += ix0 + ((ix1&sign)>>31);
   1.142 +	    ix1 += ix1;
   1.143 +	    r>>=1;
   1.144 +	}
   1.145  
   1.146 -    r = sign;
   1.147 -    while (r != 0) {
   1.148 -        t1 = s1 + r;
   1.149 -        t = s0;
   1.150 -        if ((t < ix0) || ((t == ix0) && (t1 <= ix1))) {
   1.151 -            s1 = t1 + r;
   1.152 -            if (((t1 & sign) == sign) && (s1 & sign) == 0)
   1.153 -                s0 += 1;
   1.154 -            ix0 -= t;
   1.155 -            if (ix1 < t1)
   1.156 -                ix0 -= 1;
   1.157 -            ix1 -= t1;
   1.158 -            q1 += r;
   1.159 -        }
   1.160 -        ix0 += ix0 + ((ix1 & sign) >> 31);
   1.161 -        ix1 += ix1;
   1.162 -        r >>= 1;
   1.163 -    }
   1.164 +	r = sign;
   1.165 +	while(r!=0) {
   1.166 +	    t1 = s1+r;
   1.167 +	    t  = s0;
   1.168 +	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
   1.169 +		s1  = t1+r;
   1.170 +		if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
   1.171 +		ix0 -= t;
   1.172 +		if (ix1 < t1) ix0 -= 1;
   1.173 +		ix1 -= t1;
   1.174 +		q1  += r;
   1.175 +	    }
   1.176 +	    ix0 += ix0 + ((ix1&sign)>>31);
   1.177 +	    ix1 += ix1;
   1.178 +	    r>>=1;
   1.179 +	}
   1.180  
   1.181      /* use floating add to find out rounding direction */
   1.182 -    if ((ix0 | ix1) != 0) {
   1.183 -        z = one - tiny;         /* trigger inexact flag */
   1.184 -        if (z >= one) {
   1.185 -            z = one + tiny;
   1.186 -            if (q1 == (u_int32_t) 0xffffffff) {
   1.187 -                q1 = 0;
   1.188 -                q += 1;
   1.189 -            } else if (z > one) {
   1.190 -                if (q1 == (u_int32_t) 0xfffffffe)
   1.191 -                    q += 1;
   1.192 -                q1 += 2;
   1.193 -            } else
   1.194 -                q1 += (q1 & 1);
   1.195 -        }
   1.196 -    }
   1.197 -    ix0 = (q >> 1) + 0x3fe00000;
   1.198 -    ix1 = q1 >> 1;
   1.199 -    if ((q & 1) == 1)
   1.200 -        ix1 |= sign;
   1.201 -    ix0 += (m << 20);
   1.202 -    INSERT_WORDS(z, ix0, ix1);
   1.203 -    return z;
   1.204 +	if((ix0|ix1)!=0) {
   1.205 +	    z = one-tiny; /* trigger inexact flag */
   1.206 +	    if (z>=one) {
   1.207 +	        z = one+tiny;
   1.208 +	        if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
   1.209 +		else if (z>one) {
   1.210 +		    if (q1==(u_int32_t)0xfffffffe) q+=1;
   1.211 +		    q1+=2;
   1.212 +		} else
   1.213 +	            q1 += (q1&1);
   1.214 +	    }
   1.215 +	}
   1.216 +	ix0 = (q>>1)+0x3fe00000;
   1.217 +	ix1 =  q1>>1;
   1.218 +	if ((q&1)==1) ix1 |= sign;
   1.219 +	ix0 += (m <<20);
   1.220 +	INSERT_WORDS(z,ix0,ix1);
   1.221 +	return z;
   1.222  }
   1.223  
   1.224  /*
   1.225 + * wrapper sqrt(x)
   1.226 + */
   1.227 +#ifndef _IEEE_LIBM
   1.228 +double sqrt(double x)
   1.229 +{
   1.230 +	double z = __ieee754_sqrt(x);
   1.231 +	if (_LIB_VERSION == _IEEE_ || isnan(x))
   1.232 +		return z;
   1.233 +	if (x < 0.0)
   1.234 +		return __kernel_standard(x, x, 26); /* sqrt(negative) */
   1.235 +	return z;
   1.236 +}
   1.237 +#else
   1.238 +strong_alias(__ieee754_sqrt, sqrt)
   1.239 +#endif
   1.240 +libm_hidden_def(sqrt)
   1.241 +
   1.242 +
   1.243 +/*
   1.244  Other methods  (use floating-point arithmetic)
   1.245  -------------
   1.246  (This is a copy of a drafted paper by Prof W. Kahan