From: Jakub Jelinek <jakub@redhat.com>
To: Ulrich Drepper <drepper@redhat.com>, Roland McGrath <roland@redhat.com>
Cc: Glibc hackers <libc-hacker@sources.redhat.com>
Subject: [PATCH] Fix sqrt special cases
Date: Fri, 11 Jun 2004 15:52:00 -0000 [thread overview]
Message-ID: <20040611133843.GN5191@sunsite.ms.mff.cuni.cz> (raw)
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
reply other threads:[~2004-06-11 15:52 UTC|newest]
Thread overview: [no followups] expand[flat|nested] mbox.gz Atom feed
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=20040611133843.GN5191@sunsite.ms.mff.cuni.cz \
--to=jakub@redhat.com \
--cc=drepper@redhat.com \
--cc=libc-hacker@sources.redhat.com \
--cc=roland@redhat.com \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).