Update.

2000-10-25 Ulrich Drepper <drepper@redhat.com> * sysdeps/ieee754/dbl-64/e_jn.c: Use __ieee754_sqrt instead of __sqrt. * sysdeps/ieee754/dbl-64/e_j1.c: Likewise. * sysdeps/ieee754/dbl-64/e_j0.c: Likewise. * sysdeps/ieee754/flt-32/e_j1f.c: Likewise. * sysdeps/ieee754/flt-32/e_j0f.c: Likewise.
git-mirror · Oct 25, 2000 · 1065998 · 1065998
1 parent 6a39d02
commit 1065998
Show file tree

Hide file tree

Showing 6 changed files with 67 additions and 59 deletions.
diff --git a/ChangeLog b/ChangeLog
@@ -1,3 +1,11 @@
+2000-10-25  Ulrich Drepper  <drepper@redhat.com>
+
+	* sysdeps/ieee754/dbl-64/e_jn.c: Use __ieee754_sqrt instead of __sqrt.
+	* sysdeps/ieee754/dbl-64/e_j1.c: Likewise.
+	* sysdeps/ieee754/dbl-64/e_j0.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_j1f.c: Likewise.
+	* sysdeps/ieee754/flt-32/e_j0f.c: Likewise.
+
 2000-10-25  David Mosberger  <davidm@hpl.hp.com>
 
 	* sysdeps/unix/sysv/linux/ia64/profil-counter.h: Multiply slot

diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c
@@ -124,10 +124,10 @@ static double zero = 0.0;
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x48000000) z = (invsqrtpi*cc)/__sqrt(x);
+		if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(x);
 		else {
 		    u = pzero(x); v = qzero(x);
-		    z = invsqrtpi*(u*cc-v*ss)/__sqrt(x);
+		    z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(x);
 		}
 		return z;
 	}
@@ -215,10 +215,10 @@ V[]  =  {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
                     if ((s*c)<zero) cc = z/ss;
                     else            ss = z/cc;
                 }
-                if(ix>0x48000000) z = (invsqrtpi*ss)/__sqrt(x);
+                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
                 else {
                     u = pzero(x); v = qzero(x);
-                    z = invsqrtpi*(u*ss+v*cc)/__sqrt(x);
+                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
                 }
                 return z;
 	}

diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c
@@ -125,10 +125,10 @@ static double zero    = 0.0;
 	 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 	 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x48000000) z = (invsqrtpi*cc)/__sqrt(y);
+		if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrt(y);
 		else {
 		    u = pone(y); v = qone(y);
-		    z = invsqrtpi*(u*cc-v*ss)/__sqrt(y);
+		    z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(y);
 		}
 		if(hx<0) return -z;
 		else  	 return  z;
@@ -214,10 +214,10 @@ static double V0[5] = {
          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
          * to compute the worse one.
          */
-                if(ix>0x48000000) z = (invsqrtpi*ss)/__sqrt(x);
+                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
                 else {
                     u = pone(x); v = qone(x);
-                    z = invsqrtpi*(u*ss+v*cc)/__sqrt(x);
+                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
                 }
                 return z;
         }

diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c
@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -18,7 +18,7 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $";
  * __ieee754_jn(n, x), __ieee754_yn(n, x)
  * floating point Bessel's function of the 1st and 2nd kind
  * of order n
- *          
+ *
  * Special cases:
  *	y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  *	y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
@@ -37,7 +37,7 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $";
  *	yn(n,x) is similar in all respects, except
  *	that forward recursion is used for all
  *	values of n>1.
- *	
+ *
  */
 
 #include "math.h"
@@ -76,7 +76,7 @@ static double zero  =  0.00000000000000000000e+00;
 	ix = 0x7fffffff&hx;
     /* if J(n,NaN) is NaN */
 	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
-	if(n<0){		
+	if(n<0){
 		n = -n;
 		x = -x;
 		hx ^= 0x80000000;
@@ -87,13 +87,13 @@ static double zero  =  0.00000000000000000000e+00;
 	x = fabs(x);
 	if((ix|lx)==0||ix>=0x7ff00000) 	/* if x is 0 or inf */
 	    b = zero;
-	else if((double)n<=x) {   
+	else if((double)n<=x) {
 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 	    if(ix>=0x52D00000) { /* x > 2**302 */
-    /* (x >> n**2) 
+    /* (x >> n**2)
      *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
      *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Let s=sin(x), c=cos(x), 
+     *	    Let s=sin(x), c=cos(x),
      *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
      *
      *		   n	sin(xn)*sqt2	cos(xn)*sqt2
@@ -109,8 +109,8 @@ static double zero  =  0.00000000000000000000e+00;
 		    case 2: temp = -__cos(x)-__sin(x); break;
 		    case 3: temp =  __cos(x)-__sin(x); break;
 		}
-		b = invsqrtpi*temp/__sqrt(x);
-	    } else {	
+		b = invsqrtpi*temp/__ieee754_sqrt(x);
+	    } else {
 	        a = __ieee754_j0(x);
 	        b = __ieee754_j1(x);
 	        for(i=1;i<n;i++){
@@ -121,7 +121,7 @@ static double zero  =  0.00000000000000000000e+00;
 	    }
 	} else {
 	    if(ix<0x3e100000) {	/* x < 2**-29 */
-    /* x is tiny, return the first Taylor expansion of J(n,x) 
+    /* x is tiny, return the first Taylor expansion of J(n,x)
      * J(n,x) = 1/n!*(x/2)^n  - ...
      */
 		if(n>33)	/* underflow */
@@ -136,14 +136,14 @@ static double zero  =  0.00000000000000000000e+00;
 		}
 	    } else {
 		/* use backward recurrence */
-		/* 			x      x^2      x^2       
+		/* 			x      x^2      x^2
 		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
 		 *			2n  - 2(n+1) - 2(n+2)
 		 *
-		 * 			1      1        1       
+		 * 			1      1        1
 		 *  (for large x)   =  ----  ------   ------   .....
 		 *			2n   2(n+1)   2(n+2)
-		 *			-- - ------ - ------ - 
+		 *			-- - ------ - ------ -
 		 *			 x     x         x
 		 *
 		 * Let w = 2n/x and h=2/x, then the above quotient
@@ -159,9 +159,9 @@ static double zero  =  0.00000000000000000000e+00;
 		 * To determine how many terms needed, let
 		 * Q(0) = w, Q(1) = w(w+h) - 1,
 		 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-		 * When Q(k) > 1e4	good for single 
-		 * When Q(k) > 1e9	good for double 
-		 * When Q(k) > 1e17	good for quadruple 
+		 * When Q(k) > 1e4	good for single
+		 * When Q(k) > 1e9	good for double
+		 * When Q(k) > 1e17	good for quadruple
 		 */
 	    /* determine k */
 		double t,v;
@@ -183,7 +183,7 @@ static double zero  =  0.00000000000000000000e+00;
 		 *  single 8.8722839355e+01
 		 *  double 7.09782712893383973096e+02
 		 *  long double 1.1356523406294143949491931077970765006170e+04
-		 *  then recurrent value may overflow and the result is 
+		 *  then recurrent value may overflow and the result is
 		 *  likely underflow to zero
 		 */
 		tmp = n;
@@ -219,9 +219,9 @@ static double zero  =  0.00000000000000000000e+00;
 }
 
 #ifdef __STDC__
-	double __ieee754_yn(int n, double x) 
+	double __ieee754_yn(int n, double x)
 #else
-	double __ieee754_yn(n,x) 
+	double __ieee754_yn(n,x)
 	int n; double x;
 #endif
 {
@@ -244,10 +244,10 @@ static double zero  =  0.00000000000000000000e+00;
 	if(n==1) return(sign*__ieee754_y1(x));
 	if(ix==0x7ff00000) return zero;
 	if(ix>=0x52D00000) { /* x > 2**302 */
-    /* (x >> n**2) 
+    /* (x >> n**2)
      *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
      *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-     *	    Let s=sin(x), c=cos(x), 
+     *	    Let s=sin(x), c=cos(x),
      *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
      *
      *		   n	sin(xn)*sqt2	cos(xn)*sqt2
@@ -263,14 +263,14 @@ static double zero  =  0.00000000000000000000e+00;
 		    case 2: temp = -__sin(x)+__cos(x); break;
 		    case 3: temp =  __sin(x)+__cos(x); break;
 		}
-		b = invsqrtpi*temp/__sqrt(x);
+		b = invsqrtpi*temp/__ieee754_sqrt(x);
 	} else {
 	    u_int32_t high;
 	    a = __ieee754_y0(x);
 	    b = __ieee754_y1(x);
 	/* quit if b is -inf */
 	    GET_HIGH_WORD(high,b);
-	    for(i=1;i<n&&high!=0xfff00000;i++){ 
+	    for(i=1;i<n&&high!=0xfff00000;i++){
 		temp = b;
 		b = ((double)(i+i)/x)*b - a;
 		GET_HIGH_WORD(high,b);

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -8,7 +8,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -27,9 +27,9 @@ static float pzerof(), qzerof();
 #endif
 
 #ifdef __STDC__
-static const float 
+static const float
 #else
-static float 
+static float
 #endif
 huge 	= 1e30,
 one	= 1.0,
@@ -52,9 +52,9 @@ static float zero = 0.0;
 #endif
 
 #ifdef __STDC__
-	float __ieee754_j0f(float x) 
+	float __ieee754_j0f(float x)
 #else
-	float __ieee754_j0f(x) 
+	float __ieee754_j0f(x)
 	float x;
 #endif
 {
@@ -79,10 +79,10 @@ static float zero = 0.0;
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x48000000) z = (invsqrtpi*cc)/__sqrtf(x);
+		if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);
 		else {
 		    u = pzerof(x); v = qzerof(x);
-		    z = invsqrtpi*(u*cc-v*ss)/__sqrtf(x);
+		    z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(x);
 		}
 		return z;
 	}
@@ -121,9 +121,9 @@ v03  =  2.5915085189e-07, /* 0x348b216c */
 v04  =  4.4111031494e-10; /* 0x2ff280c2 */
 
 #ifdef __STDC__
-	float __ieee754_y0f(float x) 
+	float __ieee754_y0f(float x)
 #else
-	float __ieee754_y0f(x) 
+	float __ieee754_y0f(x)
 	float x;
 #endif
 {
@@ -133,7 +133,7 @@ v04  =  4.4111031494e-10; /* 0x2ff280c2 */
 	GET_FLOAT_WORD(hx,x);
         ix = 0x7fffffff&hx;
     /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7f800000) return  one/(x+x*x); 
+	if(ix>=0x7f800000) return  one/(x+x*x);
         if(ix==0) return -one/zero;
         if(hx<0) return zero/zero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
@@ -161,10 +161,10 @@ v04  =  4.4111031494e-10; /* 0x2ff280c2 */
                     if ((s*c)<zero) cc = z/ss;
                     else            ss = z/cc;
                 }
-                if(ix>0x48000000) z = (invsqrtpi*ss)/__sqrtf(x);
+                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
                 else {
                     u = pzerof(x); v = qzerof(x);
-                    z = invsqrtpi*(u*ss+v*cc)/__sqrtf(x);
+                    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x);
                 }
                 return z;
 	}
@@ -306,7 +306,7 @@ static float pS2[5] = {
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
 	return one+ r/s;
 }
-		
+
 
 /* For x >= 8, the asymptotic expansions of qzero is
  *	-1/8 s + 75/1024 s^3 - ..., where s = 1/x.