Index: head/lib/msun/src/e_j0f.c
===================================================================
--- head/lib/msun/src/e_j0f.c	(revision 279490)
+++ head/lib/msun/src/e_j0f.c	(revision 279491)
@@ -1,337 +1,337 @@
 /* e_j0f.c -- float version of e_j0.c.
  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
  */
 
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
  * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
 
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
 #include "math.h"
 #include "math_private.h"
 
 static float pzerof(float), qzerof(float);
 
 static const float
 huge 	= 1e30,
 one	= 1.0,
 invsqrtpi=  5.6418961287e-01, /* 0x3f106ebb */
 tpi      =  6.3661974669e-01, /* 0x3f22f983 */
  		/* R0/S0 on [0, 2.00] */
 R02  =  1.5625000000e-02, /* 0x3c800000 */
 R03  = -1.8997929874e-04, /* 0xb947352e */
 R04  =  1.8295404516e-06, /* 0x35f58e88 */
 R05  = -4.6183270541e-09, /* 0xb19eaf3c */
 S01  =  1.5619102865e-02, /* 0x3c7fe744 */
 S02  =  1.1692678527e-04, /* 0x38f53697 */
 S03  =  5.1354652442e-07, /* 0x3509daa6 */
 S04  =  1.1661400734e-09; /* 0x30a045e8 */
 
 static const float zero = 0.0;
 
 float
 __ieee754_j0f(float x)
 {
 	float z, s,c,ss,cc,r,u,v;
 	int32_t hx,ix;
 
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
 	if(ix>=0x7f800000) return one/(x*x);
 	x = fabsf(x);
 	if(ix >= 0x40000000) {	/* |x| >= 2.0 */
 		s = sinf(x);
 		c = cosf(x);
 		ss = s-c;
 		cc = s+c;
 		if(ix<0x7f000000) {  /* make sure x+x not overflow */
 		    z = -cosf(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>0x54000000) z = (invsqrtpi*cc)/sqrtf(x);
+		if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */
 		else {
 		    u = pzerof(x); v = qzerof(x);
 		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
 		}
 		return z;
 	}
 	if(ix<0x3b000000) {	/* |x| < 2**-9 */
 	    if(huge+x>one) {	/* raise inexact if x != 0 */
 	        if(ix<0x39800000) return one;	/* |x|<2**-12 */
 	        else 	      return one - x*x/4;
 	    }
 	}
 	z = x*x;
 	r =  z*(R02+z*(R03+z*(R04+z*R05)));
 	s =  one+z*(S01+z*(S02+z*(S03+z*S04)));
 	if(ix < 0x3F800000) {	/* |x| < 1.00 */
 	    return one + z*((float)-0.25+(r/s));
 	} else {
 	    u = (float)0.5*x;
 	    return((one+u)*(one-u)+z*(r/s));
 	}
 }
 
 static const float
 u00  = -7.3804296553e-02, /* 0xbd9726b5 */
 u01  =  1.7666645348e-01, /* 0x3e34e80d */
 u02  = -1.3818567619e-02, /* 0xbc626746 */
 u03  =  3.4745343146e-04, /* 0x39b62a69 */
 u04  = -3.8140706238e-06, /* 0xb67ff53c */
 u05  =  1.9559013964e-08, /* 0x32a802ba */
 u06  = -3.9820518410e-11, /* 0xae2f21eb */
 v01  =  1.2730483897e-02, /* 0x3c509385 */
 v02  =  7.6006865129e-05, /* 0x389f65e0 */
 v03  =  2.5915085189e-07, /* 0x348b216c */
 v04  =  4.4111031494e-10; /* 0x2ff280c2 */
 
 float
 __ieee754_y0f(float x)
 {
 	float z, s,c,ss,cc,u,v;
 	int32_t hx,ix;
 
 	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==0) return -one/zero;
         if(hx<0) return zero/zero;
         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.
          */
                 s = sinf(x);
                 c = cosf(x);
                 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<0x7f000000) {  /* make sure x+x not overflow */
                     z = -cosf(x+x);
                     if ((s*c)<zero) cc = z/ss;
                     else            ss = z/cc;
                 }
-                if(ix>0x54800000) z = (invsqrtpi*ss)/sqrtf(x);
+                if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
                 else {
                     u = pzerof(x); v = qzerof(x);
                     z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
                 }
                 return z;
 	}
-	if(ix<=0x32000000) {	/* x < 2**-27 */
+	if(ix<=0x39000000) {	/* x < 2**-13 */
 	    return(u00 + tpi*__ieee754_logf(x));
 	}
 	z = x*x;
 	u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
 	v = one+z*(v01+z*(v02+z*(v03+z*v04)));
 	return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x)));
 }
 
 /* The asymptotic expansions of pzero is
  *	1 - 9/128 s^2 + 11025/98304 s^4 - ...,	where s = 1/x.
  * For x >= 2, We approximate pzero by
  * 	pzero(x) = 1 + (R/S)
  * where  R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
  * 	  S = 1 + pS0*s^2 + ... + pS4*s^10
  * and
  *	| pzero(x)-1-R/S | <= 2  ** ( -60.26)
  */
 static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
   0.0000000000e+00, /* 0x00000000 */
  -7.0312500000e-02, /* 0xbd900000 */
  -8.0816707611e+00, /* 0xc1014e86 */
  -2.5706311035e+02, /* 0xc3808814 */
  -2.4852163086e+03, /* 0xc51b5376 */
  -5.2530439453e+03, /* 0xc5a4285a */
 };
 static const float pS8[5] = {
   1.1653436279e+02, /* 0x42e91198 */
   3.8337448730e+03, /* 0x456f9beb */
   4.0597855469e+04, /* 0x471e95db */
   1.1675296875e+05, /* 0x47e4087c */
   4.7627726562e+04, /* 0x473a0bba */
 };
 static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  -1.1412546255e-11, /* 0xad48c58a */
  -7.0312492549e-02, /* 0xbd8fffff */
  -4.1596107483e+00, /* 0xc0851b88 */
  -6.7674766541e+01, /* 0xc287597b */
  -3.3123129272e+02, /* 0xc3a59d9b */
  -3.4643338013e+02, /* 0xc3ad3779 */
 };
 static const float pS5[5] = {
   6.0753936768e+01, /* 0x42730408 */
   1.0512523193e+03, /* 0x44836813 */
   5.9789707031e+03, /* 0x45bad7c4 */
   9.6254453125e+03, /* 0x461665c8 */
   2.4060581055e+03, /* 0x451660ee */
 };
 
 static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
  -2.5470459075e-09, /* 0xb12f081b */
  -7.0311963558e-02, /* 0xbd8fffb8 */
  -2.4090321064e+00, /* 0xc01a2d95 */
  -2.1965976715e+01, /* 0xc1afba52 */
  -5.8079170227e+01, /* 0xc2685112 */
  -3.1447946548e+01, /* 0xc1fb9565 */
 };
 static const float pS3[5] = {
   3.5856033325e+01, /* 0x420f6c94 */
   3.6151397705e+02, /* 0x43b4c1ca */
   1.1936077881e+03, /* 0x44953373 */
   1.1279968262e+03, /* 0x448cffe6 */
   1.7358093262e+02, /* 0x432d94b8 */
 };
 
 static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  -8.8753431271e-08, /* 0xb3be98b7 */
  -7.0303097367e-02, /* 0xbd8ffb12 */
  -1.4507384300e+00, /* 0xbfb9b1cc */
  -7.6356959343e+00, /* 0xc0f4579f */
  -1.1193166733e+01, /* 0xc1331736 */
  -3.2336456776e+00, /* 0xc04ef40d */
 };
 static const float pS2[5] = {
   2.2220300674e+01, /* 0x41b1c32d */
   1.3620678711e+02, /* 0x430834f0 */
   2.7047027588e+02, /* 0x43873c32 */
   1.5387539673e+02, /* 0x4319e01a */
   1.4657617569e+01, /* 0x416a859a */
 };
 
 	static float pzerof(float x)
 {
 	const float *p,*q;
 	float z,r,s;
 	int32_t ix;
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
 	if(ix>=0x41000000)     {p = pR8; q= pS8;}
-	else if(ix>=0x40f71c58){p = pR5; q= pS5;}
-	else if(ix>=0x4036db68){p = pR3; q= pS3;}
+	else if(ix>=0x409173eb){p = pR5; q= pS5;}
+	else if(ix>=0x4036d917){p = pR3; q= pS3;}
 	else                   {p = pR2; q= pS2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[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.
  * We approximate pzero by
  * 	qzero(x) = s*(-1.25 + (R/S))
  * where  R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
  * 	  S = 1 + qS0*s^2 + ... + qS5*s^12
  * and
  *	| qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
  */
 static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
   0.0000000000e+00, /* 0x00000000 */
   7.3242187500e-02, /* 0x3d960000 */
   1.1768206596e+01, /* 0x413c4a93 */
   5.5767340088e+02, /* 0x440b6b19 */
   8.8591972656e+03, /* 0x460a6cca */
   3.7014625000e+04, /* 0x471096a0 */
 };
 static const float qS8[6] = {
   1.6377603149e+02, /* 0x4323c6aa */
   8.0983447266e+03, /* 0x45fd12c2 */
   1.4253829688e+05, /* 0x480b3293 */
   8.0330925000e+05, /* 0x49441ed4 */
   8.4050156250e+05, /* 0x494d3359 */
  -3.4389928125e+05, /* 0xc8a7eb69 */
 };
 
 static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
   1.8408595828e-11, /* 0x2da1ec79 */
   7.3242180049e-02, /* 0x3d95ffff */
   5.8356351852e+00, /* 0x40babd86 */
   1.3511157227e+02, /* 0x43071c90 */
   1.0272437744e+03, /* 0x448067cd */
   1.9899779053e+03, /* 0x44f8bf4b */
 };
 static const float qS5[6] = {
   8.2776611328e+01, /* 0x42a58da0 */
   2.0778142090e+03, /* 0x4501dd07 */
   1.8847289062e+04, /* 0x46933e94 */
   5.6751113281e+04, /* 0x475daf1d */
   3.5976753906e+04, /* 0x470c88c1 */
  -5.3543427734e+03, /* 0xc5a752be */
 };
 
 static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
   4.3774099900e-09, /* 0x3196681b */
   7.3241114616e-02, /* 0x3d95ff70 */
   3.3442313671e+00, /* 0x405607e3 */
   4.2621845245e+01, /* 0x422a7cc5 */
   1.7080809021e+02, /* 0x432acedf */
   1.6673394775e+02, /* 0x4326bbe4 */
 };
 static const float qS3[6] = {
   4.8758872986e+01, /* 0x42430916 */
   7.0968920898e+02, /* 0x44316c1c */
   3.7041481934e+03, /* 0x4567825f */
   6.4604252930e+03, /* 0x45c9e367 */
   2.5163337402e+03, /* 0x451d4557 */
  -1.4924745178e+02, /* 0xc3153f59 */
 };
 
 static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
   1.5044444979e-07, /* 0x342189db */
   7.3223426938e-02, /* 0x3d95f62a */
   1.9981917143e+00, /* 0x3fffc4bf */
   1.4495602608e+01, /* 0x4167edfd */
   3.1666231155e+01, /* 0x41fd5471 */
   1.6252708435e+01, /* 0x4182058c */
 };
 static const float qS2[6] = {
   3.0365585327e+01, /* 0x41f2ecb8 */
   2.6934811401e+02, /* 0x4386ac8f */
   8.4478375244e+02, /* 0x44533229 */
   8.8293585205e+02, /* 0x445cbbe5 */
   2.1266638184e+02, /* 0x4354aa98 */
  -5.3109550476e+00, /* 0xc0a9f358 */
 };
 
 	static float qzerof(float x)
 {
 	const float *p,*q;
 	float s,r,z;
 	int32_t ix;
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
 	if(ix>=0x41000000)     {p = qR8; q= qS8;}
-	else if(ix>=0x40f71c58){p = qR5; q= qS5;}
-	else if(ix>=0x4036db68){p = qR3; q= qS3;}
+	else if(ix>=0x409173eb){p = qR5; q= qS5;}
+	else if(ix>=0x4036d917){p = qR3; q= qS3;}
 	else                   {p = qR2; q= qS2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
 	return (-(float).125 + r/s)/x;
 }
Index: head/lib/msun/src/e_j1f.c
===================================================================
--- head/lib/msun/src/e_j1f.c	(revision 279490)
+++ head/lib/msun/src/e_j1f.c	(revision 279491)
@@ -1,333 +1,333 @@
 /* e_j1f.c -- float version of e_j1.c.
  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
  */
 
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
  * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
 
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
 #include "math.h"
 #include "math_private.h"
 
 static float ponef(float), qonef(float);
 
 static const float
 huge    = 1e30,
 one	= 1.0,
 invsqrtpi=  5.6418961287e-01, /* 0x3f106ebb */
 tpi      =  6.3661974669e-01, /* 0x3f22f983 */
 	/* R0/S0 on [0,2] */
 r00  = -6.2500000000e-02, /* 0xbd800000 */
 r01  =  1.4070566976e-03, /* 0x3ab86cfd */
 r02  = -1.5995563444e-05, /* 0xb7862e36 */
 r03  =  4.9672799207e-08, /* 0x335557d2 */
 s01  =  1.9153760746e-02, /* 0x3c9ce859 */
 s02  =  1.8594678841e-04, /* 0x3942fab6 */
 s03  =  1.1771846857e-06, /* 0x359dffc2 */
 s04  =  5.0463624390e-09, /* 0x31ad6446 */
 s05  =  1.2354227016e-11; /* 0x2d59567e */
 
 static const float zero    = 0.0;
 
 float
 __ieee754_j1f(float x)
 {
 	float z, s,c,ss,cc,r,u,v,y;
 	int32_t hx,ix;
 
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
 	if(ix>=0x7f800000) return one/x;
 	y = fabsf(x);
 	if(ix >= 0x40000000) {	/* |x| >= 2.0 */
 		s = sinf(y);
 		c = cosf(y);
 		ss = -s-c;
 		cc = s-c;
 		if(ix<0x7f000000) {  /* make sure y+y not overflow */
 		    z = cosf(y+y);
 		    if ((s*c)>zero) cc = z/ss;
 		    else 	    ss = z/cc;
 		}
 	/*
 	 * 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>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
+		if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */
 		else {
 		    u = ponef(y); v = qonef(y);
 		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
 		}
 		if(hx<0) return -z;
 		else  	 return  z;
 	}
-	if(ix<0x32000000) {	/* |x|<2**-27 */
+	if(ix<0x39000000) {	/* |x|<2**-13 */
 	    if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */
 	}
 	z = x*x;
 	r =  z*(r00+z*(r01+z*(r02+z*r03)));
 	s =  one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
 	r *= x;
 	return(x*(float)0.5+r/s);
 }
 
 static const float U0[5] = {
  -1.9605709612e-01, /* 0xbe48c331 */
   5.0443872809e-02, /* 0x3d4e9e3c */
  -1.9125689287e-03, /* 0xbafaaf2a */
   2.3525259166e-05, /* 0x37c5581c */
  -9.1909917899e-08, /* 0xb3c56003 */
 };
 static const float V0[5] = {
   1.9916731864e-02, /* 0x3ca3286a */
   2.0255257550e-04, /* 0x3954644b */
   1.3560879779e-06, /* 0x35b602d4 */
   6.2274145840e-09, /* 0x31d5f8eb */
   1.6655924903e-11, /* 0x2d9281cf */
 };
 
 float
 __ieee754_y1f(float x)
 {
 	float z, s,c,ss,cc,u,v;
 	int32_t hx,ix;
 
 	GET_FLOAT_WORD(hx,x);
         ix = 0x7fffffff&hx;
     /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 	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 */
                 s = sinf(x);
                 c = cosf(x);
                 ss = -s-c;
                 cc = s-c;
                 if(ix<0x7f000000) {  /* make sure x+x not overflow */
                     z = cosf(x+x);
                     if ((s*c)>zero) cc = z/ss;
                     else            ss = z/cc;
                 }
         /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
          * where x0 = x-3pi/4
          *      Better formula:
          *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
          *                      =  1/sqrt(2) * (sin(x) - cos(x))
          *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
          *                      = -1/sqrt(2) * (cos(x) + sin(x))
          * To avoid cancellation, use
          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
          * to compute the worse one.
          */
-                if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
+                if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
                 else {
                     u = ponef(x); v = qonef(x);
                     z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
                 }
                 return z;
         }
-        if(ix<=0x24800000) {    /* x < 2**-54 */
+        if(ix<=0x33000000) {    /* x < 2**-25 */
             return(-tpi/x);
         }
         z = x*x;
         u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
         v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
         return(x*(u/v) + tpi*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x));
 }
 
 /* For x >= 8, the asymptotic expansions of pone is
  *	1 + 15/128 s^2 - 4725/2^15 s^4 - ...,	where s = 1/x.
  * We approximate pone by
  * 	pone(x) = 1 + (R/S)
  * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
  * 	  S = 1 + ps0*s^2 + ... + ps4*s^10
  * and
  *	| pone(x)-1-R/S | <= 2  ** ( -60.06)
  */
 
 static const float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
   0.0000000000e+00, /* 0x00000000 */
   1.1718750000e-01, /* 0x3df00000 */
   1.3239480972e+01, /* 0x4153d4ea */
   4.1205184937e+02, /* 0x43ce06a3 */
   3.8747453613e+03, /* 0x45722bed */
   7.9144794922e+03, /* 0x45f753d6 */
 };
 static const float ps8[5] = {
   1.1420736694e+02, /* 0x42e46a2c */
   3.6509309082e+03, /* 0x45642ee5 */
   3.6956207031e+04, /* 0x47105c35 */
   9.7602796875e+04, /* 0x47bea166 */
   3.0804271484e+04, /* 0x46f0a88b */
 };
 
 static const float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
   1.3199052094e-11, /* 0x2d68333f */
   1.1718749255e-01, /* 0x3defffff */
   6.8027510643e+00, /* 0x40d9b023 */
   1.0830818176e+02, /* 0x42d89dca */
   5.1763616943e+02, /* 0x440168b7 */
   5.2871520996e+02, /* 0x44042dc6 */
 };
 static const float ps5[5] = {
   5.9280597687e+01, /* 0x426d1f55 */
   9.9140142822e+02, /* 0x4477d9b1 */
   5.3532670898e+03, /* 0x45a74a23 */
   7.8446904297e+03, /* 0x45f52586 */
   1.5040468750e+03, /* 0x44bc0180 */
 };
 
 static const float pr3[6] = {
   3.0250391081e-09, /* 0x314fe10d */
   1.1718686670e-01, /* 0x3defffab */
   3.9329774380e+00, /* 0x407bb5e7 */
   3.5119403839e+01, /* 0x420c7a45 */
   9.1055007935e+01, /* 0x42b61c2a */
   4.8559066772e+01, /* 0x42423c7c */
 };
 static const float ps3[5] = {
   3.4791309357e+01, /* 0x420b2a4d */
   3.3676245117e+02, /* 0x43a86198 */
   1.0468714600e+03, /* 0x4482dbe3 */
   8.9081134033e+02, /* 0x445eb3ed */
   1.0378793335e+02, /* 0x42cf936c */
 };
 
 static const float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
   1.0771083225e-07, /* 0x33e74ea8 */
   1.1717621982e-01, /* 0x3deffa16 */
   2.3685150146e+00, /* 0x401795c0 */
   1.2242610931e+01, /* 0x4143e1bc */
   1.7693971634e+01, /* 0x418d8d41 */
   5.0735230446e+00, /* 0x40a25a4d */
 };
 static const float ps2[5] = {
   2.1436485291e+01, /* 0x41ab7dec */
   1.2529022980e+02, /* 0x42fa9499 */
   2.3227647400e+02, /* 0x436846c7 */
   1.1767937469e+02, /* 0x42eb5bd7 */
   8.3646392822e+00, /* 0x4105d590 */
 };
 
 	static float ponef(float x)
 {
 	const float *p,*q;
 	float z,r,s;
         int32_t ix;
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
         if(ix>=0x41000000)     {p = pr8; q= ps8;}
-        else if(ix>=0x40f71c58){p = pr5; q= ps5;}
-        else if(ix>=0x4036db68){p = pr3; q= ps3;}
+        else if(ix>=0x409173eb){p = pr5; q= ps5;}
+        else if(ix>=0x4036d917){p = pr3; q= ps3;}
 	else                   {p = pr2; q= ps2;}	/* ix>=0x40000000 */
         z = one/(x*x);
         r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[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 qone is
  *	3/8 s - 105/1024 s^3 - ..., where s = 1/x.
  * We approximate pone by
  * 	qone(x) = s*(0.375 + (R/S))
  * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
  * 	  S = 1 + qs1*s^2 + ... + qs6*s^12
  * and
  *	| qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
  */
 
 static const float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
   0.0000000000e+00, /* 0x00000000 */
  -1.0253906250e-01, /* 0xbdd20000 */
  -1.6271753311e+01, /* 0xc1822c8d */
  -7.5960174561e+02, /* 0xc43de683 */
  -1.1849806641e+04, /* 0xc639273a */
  -4.8438511719e+04, /* 0xc73d3683 */
 };
 static const float qs8[6] = {
   1.6139537048e+02, /* 0x43216537 */
   7.8253862305e+03, /* 0x45f48b17 */
   1.3387534375e+05, /* 0x4802bcd6 */
   7.1965775000e+05, /* 0x492fb29c */
   6.6660125000e+05, /* 0x4922be94 */
  -2.9449025000e+05, /* 0xc88fcb48 */
 };
 
 static const float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
  -2.0897993405e-11, /* 0xadb7d219 */
  -1.0253904760e-01, /* 0xbdd1fffe */
  -8.0564479828e+00, /* 0xc100e736 */
  -1.8366960144e+02, /* 0xc337ab6b */
  -1.3731937256e+03, /* 0xc4aba633 */
  -2.6124443359e+03, /* 0xc523471c */
 };
 static const float qs5[6] = {
   8.1276550293e+01, /* 0x42a28d98 */
   1.9917987061e+03, /* 0x44f8f98f */
   1.7468484375e+04, /* 0x468878f8 */
   4.9851425781e+04, /* 0x4742bb6d */
   2.7948074219e+04, /* 0x46da5826 */
  -4.7191835938e+03, /* 0xc5937978 */
 };
 
 static const float qr3[6] = {
  -5.0783124372e-09, /* 0xb1ae7d4f */
  -1.0253783315e-01, /* 0xbdd1ff5b */
  -4.6101160049e+00, /* 0xc0938612 */
  -5.7847221375e+01, /* 0xc267638e */
  -2.2824453735e+02, /* 0xc3643e9a */
  -2.1921012878e+02, /* 0xc35b35cb */
 };
 static const float qs3[6] = {
   4.7665153503e+01, /* 0x423ea91e */
   6.7386511230e+02, /* 0x4428775e */
   3.3801528320e+03, /* 0x45534272 */
   5.5477290039e+03, /* 0x45ad5dd5 */
   1.9031191406e+03, /* 0x44ede3d0 */
  -1.3520118713e+02, /* 0xc3073381 */
 };
 
 static const float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
  -1.7838172539e-07, /* 0xb43f8932 */
  -1.0251704603e-01, /* 0xbdd1f475 */
  -2.7522056103e+00, /* 0xc0302423 */
  -1.9663616180e+01, /* 0xc19d4f16 */
  -4.2325313568e+01, /* 0xc2294d1f */
  -2.1371921539e+01, /* 0xc1aaf9b2 */
 };
 static const float qs2[6] = {
   2.9533363342e+01, /* 0x41ec4454 */
   2.5298155212e+02, /* 0x437cfb47 */
   7.5750280762e+02, /* 0x443d602e */
   7.3939318848e+02, /* 0x4438d92a */
   1.5594900513e+02, /* 0x431bf2f2 */
  -4.9594988823e+00, /* 0xc09eb437 */
 };
 
 	static float qonef(float x)
 {
 	const float *p,*q;
 	float  s,r,z;
 	int32_t ix;
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qr8; q= qs8;}
-	else if(ix>=0x40f71c58){p = qr5; q= qs5;}
-	else if(ix>=0x4036db68){p = qr3; q= qs3;}
+	if(ix>=0x41000000)     {p = qr8; q= qs8;}
+	else if(ix>=0x409173eb){p = qr5; q= qs5;}
+	else if(ix>=0x4036d917){p = qr3; q= qs3;}
 	else                   {p = qr2; q= qs2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
 	return ((float).375 + r/s)/x;
 }