diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_sqrt.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_sqrt.c | 95 |
1 files changed, 52 insertions, 43 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c index 54610eeeab..854ae38c41 100644 --- a/sysdeps/ieee754/dbl-64/e_sqrt.c +++ b/sysdeps/ieee754/dbl-64/e_sqrt.c @@ -44,58 +44,67 @@ /* it computes the correctly rounded (to nearest) value of square */ /* root of x. */ /*********************************************************************/ -double __ieee754_sqrt(double x) { +double +__ieee754_sqrt (double x) +{ #include "uroot.h" static const double rt0 = 9.99999999859990725855365213134618E-01, rt1 = 4.99999999495955425917856814202739E-01, rt2 = 3.75017500867345182581453026130850E-01, rt3 = 3.12523626554518656309172508769531E-01; - static const double big = 134217728.0; - double y,t,del,res,res1,hy,z,zz,p,hx,tx,ty,s; - mynumber a,c={{0,0}}; + static const double big = 134217728.0; + double y, t, del, res, res1, hy, z, zz, p, hx, tx, ty, s; + mynumber a, c = { { 0, 0 } }; int4 k; - a.x=x; - k=a.i[HIGH_HALF]; - a.i[HIGH_HALF]=(k&0x001fffff)|0x3fe00000; - t=inroot[(k&0x001fffff)>>14]; - s=a.x; + a.x = x; + k = a.i[HIGH_HALF]; + a.i[HIGH_HALF] = (k & 0x001fffff) | 0x3fe00000; + t = inroot[(k & 0x001fffff) >> 14]; + s = a.x; /*----------------- 2^-1022 <= | x |< 2^1024 -----------------*/ - if (k>0x000fffff && k<0x7ff00000) { - fenv_t env; - libc_feholdexcept (&env); - double ret; - y=1.0-t*(t*s); - t=t*(rt0+y*(rt1+y*(rt2+y*rt3))); - c.i[HIGH_HALF]=0x20000000+((k&0x7fe00000)>>1); - y=t*s; - hy=(y+big)-big; - del=0.5*t*((s-hy*hy)-(y-hy)*(y+hy)); - res=y+del; - if (res == (res+1.002*((y-res)+del))) ret = res*c.x; - else { - res1=res+1.5*((y-res)+del); - EMULV(res,res1,z,zz,p,hx,tx,hy,ty); /* (z+zz)=res*res1 */ - ret = ((((z-s)+zz)<0)?max(res,res1):min(res,res1))*c.x; + if (k > 0x000fffff && k < 0x7ff00000) + { + fenv_t env; + libc_feholdexcept (&env); + double ret; + y = 1.0 - t * (t * s); + t = t * (rt0 + y * (rt1 + y * (rt2 + y * rt3))); + c.i[HIGH_HALF] = 0x20000000 + ((k & 0x7fe00000) >> 1); + y = t * s; + hy = (y + big) - big; + del = 0.5 * t * ((s - hy * hy) - (y - hy) * (y + hy)); + res = y + del; + if (res == (res + 1.002 * ((y - res) + del))) + ret = res * c.x; + else + { + res1 = res + 1.5 * ((y - res) + del); + EMULV (res, res1, z, zz, p, hx, tx, hy, ty); /* (z+zz)=res*res1 */ + ret = ((((z - s) + zz) < 0) ? max (res, res1) : + min (res, res1)) * c.x; + } + math_force_eval (ret); + libc_fesetenv (&env); + if (x / ret != ret) + { + double force_inexact = 1.0 / 3.0; + math_force_eval (force_inexact); + } + /* Otherwise (x / ret == ret), either the square root was exact or + the division was inexact. */ + return ret; + } + else + { + if ((k & 0x7ff00000) == 0x7ff00000) + return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ + if (x == 0) + return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */ + if (k < 0) + return (x - x) / (x - x); /* sqrt(-ve)=sNaN */ + return tm256.x * __ieee754_sqrt (x * t512.x); } - math_force_eval (ret); - libc_fesetenv (&env); - if (x / ret != ret) - { - double force_inexact = 1.0 / 3.0; - math_force_eval (force_inexact); - } - /* Otherwise (x / ret == ret), either the square root was exact or - the division was inexact. */ - return ret; - } - else { - if ((k & 0x7ff00000) == 0x7ff00000) - return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ - if (x==0) return x; /* sqrt(+0)=+0, sqrt(-0)=-0 */ - if (k<0) return (x-x)/(x-x); /* sqrt(-ve)=sNaN */ - return tm256.x*__ieee754_sqrt(x*t512.x); - } } strong_alias (__ieee754_sqrt, __sqrt_finite) |