aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorOndřej Bílka <neleai@seznam.cz>2013-10-17 16:03:24 +0200
committerOndřej Bílka <neleai@seznam.cz>2013-10-17 16:03:24 +0200
commitc5d5d574cbfa96d0f6c1db24d1e072c472627e41 (patch)
tree83b97e29ee65636dfe1247ea8d2344ca3f0b04b4 /sysdeps/ieee754/dbl-64/e_sqrt.c
parentsoft-fp: make extensions quiet signaling NaNs (bug 16041). (diff)
downloadglibc-c5d5d574cbfa96d0f6c1db24d1e072c472627e41.tar.gz
glibc-c5d5d574cbfa96d0f6c1db24d1e072c472627e41.tar.bz2
glibc-c5d5d574cbfa96d0f6c1db24d1e072c472627e41.zip
Format floating routines.
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_sqrt.c')
-rw-r--r--sysdeps/ieee754/dbl-64/e_sqrt.c95
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)