This is the mail archive of the libc-hacker@sources.redhat.com mailing list for the glibc project.
Note that libc-hacker is a closed list. You may look at the archives of this list, but subscription and posting are not open.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Other format: | [Raw text] |
Hi! GCC 3.4+ (unlike 3.3 and 3.2) will not optimize the division in: double foo (void) { static const double big = 134217728.0, big1 = 134217729.0; return (big1-big1)/(big-big); } out (IMHO 3.4+ is correct here, unless -ffast-math the division needs to be done, so that the exception is generated at runtime). But e_sqrt after IBM rewrite is apparently relying on this being optimized out (even in that case it is not correct, since sqrt(-inf) or sqrt(-ve) should generate exception and a non-canonical NaN can have negative sign bit set). This patch puts in what Sun e_sqrt was doing for the special cases. 2004-06-11 Jakub Jelinek <jakub@redhat.com> * sysdeps/ieee754/dbl-64/e_sqrt.c (__ieee754_sqrt): Handle special cases properly. --- libc/sysdeps/ieee754/dbl-64/e_sqrt.c.jj 2002-08-26 18:40:37.000000000 -0400 +++ libc/sysdeps/ieee754/dbl-64/e_sqrt.c 2004-06-11 11:17:02.000000000 -0400 @@ -41,7 +41,7 @@ #include "math_private.h" /*********************************************************************/ -/* An ultimate aqrt routine. Given an IEEE double machine number x */ +/* An ultimate sqrt routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of square */ /* root of x. */ /*********************************************************************/ @@ -52,7 +52,7 @@ double __ieee754_sqrt(double x) { rt1 = 4.99999999495955425917856814202739E-01, rt2 = 3.75017500867345182581453026130850E-01, rt3 = 3.12523626554518656309172508769531E-01; - static const double big = 134217728.0, big1 = 134217729.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; @@ -79,13 +79,10 @@ double __ieee754_sqrt(double x) { } } else { - if (k>0x7ff00000) /* x -> infinity */ - return (big1-big1)/(big-big); - if (k<0x00100000) { /* x -> -infinity */ - if (x==0) return x; - if (k<0) return (big1-big1)/(big-big); - else return tm256.x*__ieee754_sqrt(x*t512.x); - } - else return (a.i[LOW_HALF]==0)?x:(big1-big1)/(big-big); + 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); } } Jakub
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |