[PATCH] Fix sqrt special cases

[Jakub Jelinek <jakub at redhat dot com>]

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