Format floating routines.

1 files changed, 213 insertions, 161 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c
index ca542756fb..d165e80925 100644
--- a/sysdeps/ieee754/dbl-64/e_j0.c
+++ b/sysdeps/ieee754/dbl-64/e_j0.c
@@ -11,7 +11,7 @@
  */
 /* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/26,
    for performance improvement on pipelined processors.
-*/
+ */
 
 /* __ieee754_j0(x), __ieee754_y0(x)
  * Bessel function of the first and second kinds of order zero.
@@ -61,144 +61,166 @@
 #include <math.h>
 #include <math_private.h>
 
-static double pzero(double), qzero(double);
+static double pzero (double), qzero (double);
 
 static const double
-huge	= 1e300,
-one	= 1.0,
-invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
-tpi      =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
-		/* R0/S0 on [0, 2.00] */
-R[]  =  {0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
- -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
-  1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
- -4.61832688532103189199e-09}, /* 0xBE33D5E7, 0x73D63FCE */
-S[]  =  {0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
-  1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
-  5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
-  1.16614003333790000205e-09}; /* 0x3E1408BC, 0xF4745D8F */
+  huge = 1e300,
+  one = 1.0,
+  invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
+  tpi = 6.36619772367581382433e-01,     /* 0x3FE45F30, 0x6DC9C883 */
+/* R0/S0 on [0, 2.00] */
+  R[] = { 0.0, 0.0, 1.56249999999999947958e-02, /* 0x3F8FFFFF, 0xFFFFFFFD */
+	  -1.89979294238854721751e-04, /* 0xBF28E6A5, 0xB61AC6E9 */
+	  1.82954049532700665670e-06, /* 0x3EBEB1D1, 0x0C503919 */
+	  -4.61832688532103189199e-09 }, /* 0xBE33D5E7, 0x73D63FCE */
+  S[] = { 0.0, 1.56191029464890010492e-02, /* 0x3F8FFCE8, 0x82C8C2A4 */
+	  1.16926784663337450260e-04, /* 0x3F1EA6D2, 0xDD57DBF4 */
+	  5.13546550207318111446e-07, /* 0x3EA13B54, 0xCE84D5A9 */
+	  1.16614003333790000205e-09 }; /* 0x3E1408BC, 0xF4745D8F */
 
 static const double zero = 0.0;
 
 double
-__ieee754_j0(double x)
+__ieee754_j0 (double x)
 {
-	double z, s,c,ss,cc,r,u,v,r1,r2,s1,s2,z2,z4;
-	int32_t hx,ix;
+  double z, s, c, ss, cc, r, u, v, r1, r2, s1, s2, z2, z4;
+  int32_t hx, ix;
 
-	GET_HIGH_WORD(hx,x);
-	ix = hx&0x7fffffff;
-	if(ix>=0x7ff00000) return one/(x*x);
-	x = fabs(x);
-	if(ix >= 0x40000000) {	/* |x| >= 2.0 */
-		__sincos (x, &s, &c);
-		ss = s-c;
-		cc = s+c;
-		if(ix<0x7fe00000) {  /* make sure x+x not overflow */
-		    z = -__cos(x+x);
-		    if ((s*c)<zero) cc = z/ss;
-		    else	    ss = z/cc;
-		}
-	/*
-	 * 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)/__ieee754_sqrt(x);
-		else {
-		    u = pzero(x); v = qzero(x);
-		    z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrt(x);
-		}
-		return z;
-	}
-	if(ix<0x3f200000) {	/* |x| < 2**-13 */
-	  math_force_eval(huge+x);	/* raise inexact if x != 0 */
-	  if(ix<0x3e400000) return one;	/* |x|<2**-27 */
-	  else	      return one - 0.25*x*x;
+  GET_HIGH_WORD (hx, x);
+  ix = hx & 0x7fffffff;
+  if (ix >= 0x7ff00000)
+    return one / (x * x);
+  x = fabs (x);
+  if (ix >= 0x40000000)         /* |x| >= 2.0 */
+    {
+      __sincos (x, &s, &c);
+      ss = s - c;
+      cc = s + c;
+      if (ix < 0x7fe00000)           /* make sure x+x not overflow */
+	{
+	  z = -__cos (x + x);
+	  if ((s * c) < zero)
+	    cc = z / ss;
+	  else
+	    ss = z / cc;
 	}
-	z = x*x;
-	r1 = z*R[2]; z2=z*z;
-	r2 = R[3]+z*R[4]; z4=z2*z2;
-	r  = r1 + z2*r2 + z4*R[5];
-	s1 = one+z*S[1];
-	s2 = S[2]+z*S[3];
-	s = s1 + z2*s2 + z4*S[4];
-	if(ix < 0x3FF00000) {	/* |x| < 1.00 */
-	    return one + z*(-0.25+(r/s));
-	} else {
-	    u = 0.5*x;
-	    return((one+u)*(one-u)+z*(r/s));
+      /*
+       * 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) / __ieee754_sqrt (x);
+      else
+	{
+	  u = pzero (x); v = qzero (x);
+	  z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt (x);
 	}
+      return z;
+    }
+  if (ix < 0x3f200000)          /* |x| < 2**-13 */
+    {
+      math_force_eval (huge + x);       /* raise inexact if x != 0 */
+      if (ix < 0x3e400000)
+	return one;                     /* |x|<2**-27 */
+      else
+	return one - 0.25 * x * x;
+    }
+  z = x * x;
+  r1 = z * R[2]; z2 = z * z;
+  r2 = R[3] + z * R[4]; z4 = z2 * z2;
+  r = r1 + z2 * r2 + z4 * R[5];
+  s1 = one + z * S[1];
+  s2 = S[2] + z * S[3];
+  s = s1 + z2 * s2 + z4 * S[4];
+  if (ix < 0x3FF00000)          /* |x| < 1.00 */
+    {
+      return one + z * (-0.25 + (r / s));
+    }
+  else
+    {
+      u = 0.5 * x;
+      return ((one + u) * (one - u) + z * (r / s));
+    }
 }
 strong_alias (__ieee754_j0, __j0_finite)
 
 static const double
-U[]  = {-7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
-  1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
- -1.38185671945596898896e-02, /* 0xBF8C4CE8, 0xB16CFA97 */
-  3.47453432093683650238e-04, /* 0x3F36C54D, 0x20B29B6B */
- -3.81407053724364161125e-06, /* 0xBECFFEA7, 0x73D25CAD */
-  1.95590137035022920206e-08, /* 0x3E550057, 0x3B4EABD4 */
- -3.98205194132103398453e-11}, /* 0xBDC5E43D, 0x693FB3C8 */
-V[]  =  {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
-  7.60068627350353253702e-05, /* 0x3F13ECBB, 0xF578C6C1 */
-  2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
-  4.41110311332675467403e-10}; /* 0x3DFE5018, 0x3BD6D9EF */
+U[] = { -7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
+	 1.76666452509181115538e-01,  /* 0x3FC69D01, 0x9DE9E3FC */
+	-1.38185671945596898896e-02,  /* 0xBF8C4CE8, 0xB16CFA97 */
+	 3.47453432093683650238e-04,  /* 0x3F36C54D, 0x20B29B6B */
+	-3.81407053724364161125e-06,  /* 0xBECFFEA7, 0x73D25CAD */
+	 1.95590137035022920206e-08,  /* 0x3E550057, 0x3B4EABD4 */
+	-3.98205194132103398453e-11 }, /* 0xBDC5E43D, 0x693FB3C8 */
+V[] = { 1.27304834834123699328e-02,   /* 0x3F8A1270, 0x91C9C71A */
+	 7.60068627350353253702e-05,   /* 0x3F13ECBB, 0xF578C6C1 */
+	 2.59150851840457805467e-07,   /* 0x3E91642D, 0x7FF202FD */
+	 4.41110311332675467403e-10 }; /* 0x3DFE5018, 0x3BD6D9EF */
 
 double
-__ieee754_y0(double x)
+__ieee754_y0 (double x)
 {
-	double z, s,c,ss,cc,u,v,z2,z4,z6,u1,u2,u3,v1,v2;
-	int32_t hx,ix,lx;
+  double z, s, c, ss, cc, u, v, z2, z4, z6, u1, u2, u3, v1, v2;
+  int32_t hx, ix, lx;
 
-	EXTRACT_WORDS(hx,lx,x);
-	ix = 0x7fffffff&hx;
-    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf.  */
-	if(ix>=0x7ff00000) return  one/(x+x*x);
-	if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception.  */
-	if(hx<0) return zero/(zero*x);
-	if(ix >= 0x40000000) {  /* |x| >= 2.0 */
-	/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
-	 * where x0 = x-pi/4
-	 *      Better formula:
-	 *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
-	 *                      =  1/sqrt(2) * (sin(x) + cos(x))
-	 *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
-	 *                      =  1/sqrt(2) * (sin(x) - cos(x))
-	 * To avoid cancellation, use
-	 *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
-	 * to compute the worse one.
-	 */
-		__sincos (x, &s, &c);
-		ss = s-c;
-		cc = s+c;
-	/*
-	 * 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<0x7fe00000) {  /* make sure x+x not overflow */
-		    z = -__cos(x+x);
-		    if ((s*c)<zero) cc = z/ss;
-		    else            ss = z/cc;
-		}
-		if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrt(x);
-		else {
-		    u = pzero(x); v = qzero(x);
-		    z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrt(x);
-		}
-		return z;
+  EXTRACT_WORDS (hx, lx, x);
+  ix = 0x7fffffff & hx;
+  /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf.  */
+  if (ix >= 0x7ff00000)
+    return one / (x + x * x);
+  if ((ix | lx) == 0)
+    return -HUGE_VAL + x;                  /* -inf and overflow exception.  */
+  if (hx < 0)
+    return zero / (zero * x);
+  if (ix >= 0x40000000)         /* |x| >= 2.0 */
+    {   /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
+		 * where x0 = x-pi/4
+		 *      Better formula:
+		 *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
+		 *                      =  1/sqrt(2) * (sin(x) + cos(x))
+		 *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
+		 *                      =  1/sqrt(2) * (sin(x) - cos(x))
+		 * To avoid cancellation, use
+		 *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+		 * to compute the worse one.
+		 */
+      __sincos (x, &s, &c);
+      ss = s - c;
+      cc = s + c;
+      /*
+       * 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 < 0x7fe00000)           /* make sure x+x not overflow */
+	{
+	  z = -__cos (x + x);
+	  if ((s * c) < zero)
+	    cc = z / ss;
+	  else
+	    ss = z / cc;
 	}
-	if(ix<=0x3e400000) {	/* x < 2**-27 */
-	    return(U[0] + tpi*__ieee754_log(x));
+      if (ix > 0x48000000)
+	z = (invsqrtpi * ss) / __ieee754_sqrt (x);
+      else
+	{
+	  u = pzero (x); v = qzero (x);
+	  z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrt (x);
 	}
-	z = x*x;
-	u1 = U[0]+z*U[1]; z2=z*z;
-	u2 = U[2]+z*U[3]; z4=z2*z2;
-	u3 = U[4]+z*U[5]; z6=z4*z2;
-	u = u1 + z2*u2 + z4*u3 + z6*U[6];
-	v1 = one+z*V[0];
-	v2 = V[1]+z*V[2];
-	v = v1 + z2*v2 + z4*V[3];
-	return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
+      return z;
+    }
+  if (ix <= 0x3e400000)         /* x < 2**-27 */
+    {
+      return (U[0] + tpi * __ieee754_log (x));
+    }
+  z = x * x;
+  u1 = U[0] + z * U[1]; z2 = z * z;
+  u2 = U[2] + z * U[3]; z4 = z2 * z2;
+  u3 = U[4] + z * U[5]; z6 = z4 * z2;
+  u = u1 + z2 * u2 + z4 * u3 + z6 * U[6];
+  v1 = one + z * V[0];
+  v2 = V[1] + z * V[2];
+  v = v1 + z2 * v2 + z4 * V[3];
+  return (u / v + tpi * (__ieee754_j0 (x) * __ieee754_log (x)));
 }
 strong_alias (__ieee754_y0, __y0_finite)
 
@@ -276,28 +298,43 @@ static const double pS2[5] = {
 };
 
 static double
-pzero(double x)
+pzero (double x)
 {
-	const double *p,*q;
-	double z,r,s,z2,z4,r1,r2,r3,s1,s2,s3;
-	int32_t ix;
-	GET_HIGH_WORD(ix,x);
-	ix &= 0x7fffffff;
-	if (ix>=0x41b00000)    {return one;}
-	else if(ix>=0x40200000){p = pR8; q= pS8;}
-	else if(ix>=0x40122E8B){p = pR5; q= pS5;}
-	else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
-	else if(ix>=0x40000000){p = pR2; q= pS2;}
-	z = one/(x*x);
-	r1 = p[0]+z*p[1]; z2=z*z;
-	r2 = p[2]+z*p[3]; z4=z2*z2;
-	r3 = p[4]+z*p[5];
-	r = r1 + z2*r2 + z4*r3;
-	s1 = one+z*q[0];
-	s2 = q[1]+z*q[2];
-	s3 = q[3]+z*q[4];
-	s = s1 + z2*s2 + z4*s3;
-	return one+ r/s;
+  const double *p, *q;
+  double z, r, s, z2, z4, r1, r2, r3, s1, s2, s3;
+  int32_t ix;
+  GET_HIGH_WORD (ix, x);
+  ix &= 0x7fffffff;
+  if (ix >= 0x41b00000)
+    {
+      return one;
+    }
+  else if (ix >= 0x40200000)
+    {
+      p = pR8; q = pS8;
+    }
+  else if (ix >= 0x40122E8B)
+    {
+      p = pR5; q = pS5;
+    }
+  else if (ix >= 0x4006DB6D)
+    {
+      p = pR3; q = pS3;
+    }
+  else if (ix >= 0x40000000)
+    {
+      p = pR2; q = pS2;
+    }
+  z = one / (x * x);
+  r1 = p[0] + z * p[1]; z2 = z * z;
+  r2 = p[2] + z * p[3]; z4 = z2 * z2;
+  r3 = p[4] + z * p[5];
+  r = r1 + z2 * r2 + z4 * r3;
+  s1 = one + z * q[0];
+  s2 = q[1] + z * q[2];
+  s3 = q[3] + z * q[4];
+  s = s1 + z2 * s2 + z4 * s3;
+  return one + r / s;
 }
 
 
@@ -379,26 +416,41 @@ static const double qS2[6] = {
 };
 
 static double
-qzero(double x)
+qzero (double x)
 {
-	const double *p,*q;
-	double s,r,z,z2,z4,z6,r1,r2,r3,s1,s2,s3;
-	int32_t ix;
-	GET_HIGH_WORD(ix,x);
-	ix &= 0x7fffffff;
-	if (ix>=0x41b00000)    {return -.125/x;}
-	else if(ix>=0x40200000){p = qR8; q= qS8;}
-	else if(ix>=0x40122E8B){p = qR5; q= qS5;}
-	else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
-	else if(ix>=0x40000000){p = qR2; q= qS2;}
-	z = one/(x*x);
-	r1 = p[0]+z*p[1]; z2=z*z;
-	r2 = p[2]+z*p[3]; z4=z2*z2;
-	r3 = p[4]+z*p[5]; z6=z4*z2;
-	r= r1 + z2*r2 + z4*r3;
-	s1 = one+z*q[0];
-	s2 = q[1]+z*q[2];
-	s3 = q[3]+z*q[4];
-	s = s1 + z2*s2 + z4*s3 +z6*q[5];
-	return (-.125 + r/s)/x;
+  const double *p, *q;
+  double s, r, z, z2, z4, z6, r1, r2, r3, s1, s2, s3;
+  int32_t ix;
+  GET_HIGH_WORD (ix, x);
+  ix &= 0x7fffffff;
+  if (ix >= 0x41b00000)
+    {
+      return -.125 / x;
+    }
+  else if (ix >= 0x40200000)
+    {
+      p = qR8; q = qS8;
+    }
+  else if (ix >= 0x40122E8B)
+    {
+      p = qR5; q = qS5;
+    }
+  else if (ix >= 0x4006DB6D)
+    {
+      p = qR3; q = qS3;
+    }
+  else if (ix >= 0x40000000)
+    {
+      p = qR2; q = qS2;
+    }
+  z = one / (x * x);
+  r1 = p[0] + z * p[1]; z2 = z * z;
+  r2 = p[2] + z * p[3]; z4 = z2 * z2;
+  r3 = p[4] + z * p[5]; z6 = z4 * z2;
+  r = r1 + z2 * r2 + z4 * r3;
+  s1 = one + z * q[0];
+  s2 = q[1] + z * q[2];
+  s3 = q[3] + z * q[4];
+  s = s1 + z2 * s2 + z4 * s3 + z6 * q[5];
+  return (-.125 + r / s) / x;
 }