Remove trailing whitespace.

27 files changed, 314 insertions, 314 deletions

diff --git a/lib/libm/common/atan2.c b/lib/libm/common/atan2.c
index 958a154..b847a1d 100644
--- a/lib/libm/common/atan2.c
+++ b/lib/libm/common/atan2.c
@@ -38,21 +38,21 @@ static char sccsid[] = "@(#)atan2.c	8.1 (Berkeley) 6/4/93";
 /* ATAN2(Y,X)
  * RETURN ARG (X+iY)
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85; 
+ * CODED IN C BY K.C. NG, 1/8/85;
  * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85.
  *
  * Required system supported functions :
  *	copysign(x,y)
  *	scalb(x,y)
  *	logb(x)
- *	
+ *
  * Method :
  *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
- *	2. Reduce x to positive by (if x and y are unexceptional): 
+ *	2. Reduce x to positive by (if x and y are unexceptional):
  *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
  *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
- *	3. According to the integer k=4t+0.25 truncated , t=y/x, the argument 
- *	   is further reduced to one of the following intervals and the 
+ *	3. According to the integer k=4t+0.25 truncated , t=y/x, the argument
+ *	   is further reduced to one of the following intervals and the
  *	   arctangent of y/x is evaluated by the corresponding formula:
  *
  *         [0,7/16]	   atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
@@ -76,32 +76,32 @@ static char sccsid[] = "@(#)atan2.c	8.1 (Berkeley) 6/4/93";
  *	ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
  *
  * Accuracy:
- *	atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, 
+ *	atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded,
  *	where
  *
  *	in decimal:
- *		pi = 3.141592653589793 23846264338327 ..... 
+ *		pi = 3.141592653589793 23846264338327 .....
  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
+ *    56 bits   PI = 3.141592653589793 227020265 ..... ,
  *
  *	in hexadecimal:
  *		pi = 3.243F6A8885A308D313198A2E....
  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *	
+ *
  *	In a test run with 356,000 random argument on [-1,1] * [-1,1] on a
  *	VAX, the maximum observed error was 1.41 ulps (units of the last place)
  *	compared with (PI/pi)*(the exact ARG(x+iy)).
  *
  * Note:
  *	We use machine PI (the true pi rounded) in place of the actual
- *	value of pi for all the trig and inverse trig functions. In general, 
- *	if trig is one of sin, cos, tan, then computed trig(y) returns the 
- *	exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig 
- *	returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the 
+ *	value of pi for all the trig and inverse trig functions. In general,
+ *	if trig is one of sin, cos, tan, then computed trig(y) returns the
+ *	exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig
+ *	returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the
  *	trig functions have period PI, and trig(arctrig(x)) returns x for
  *	all critical values x.
- *	
+ *
  * Constants:
  * The hexadecimal values are the intended ones for the following constants.
  * The decimal values may be used, provided that the compiler will convert
@@ -174,7 +174,7 @@ ic(a11,    1.6438029044759730479E-2  ,  -6,  1.0D52174A1BB54)
 
 double atan2(y,x)
 double  y,x;
-{  
+{
 	static const double zero=0, one=1, small=1.0E-9, big=1.0E18;
 	double t,z,signy,signx,hi,lo;
 	int k,m;
@@ -185,8 +185,8 @@ double  y,x;
 #endif	/* !defined(vax)&&!defined(tahoe) */
 
     /* copy down the sign of y and x */
-	signy = copysign(one,y) ;  
-	signx = copysign(one,x) ;  
+	signy = copysign(one,y) ;
+	signx = copysign(one,x) ;
 
     /* if x is 1.0, goto begin */
 	if(x==1) { y=copysign(y,one); t=y; if(finite(t)) goto begin;}
@@ -196,10 +196,10 @@ double  y,x;
 
     /* when x = 0 */
 	if(x==zero) return(copysign(PIo2,signy));
-	    
+
     /* when x is INF */
 	if(!finite(x))
-	    if(!finite(y)) 
+	    if(!finite(y))
 		return(copysign((signx==one)?PIo4:3*PIo4,signy));
 	    else
 		return(copysign((signx==one)?zero:PI,signy));
@@ -208,43 +208,43 @@ double  y,x;
 	if(!finite(y)) return(copysign(PIo2,signy));
 
     /* compute y/x */
-	x=copysign(x,one); 
-	y=copysign(y,one); 
-	if((m=(k=logb(y))-logb(x)) > 60) t=big+big; 
+	x=copysign(x,one);
+	y=copysign(y,one);
+	if((m=(k=logb(y))-logb(x)) > 60) t=big+big;
 	    else if(m < -80 ) t=y/x;
 	    else { t = y/x ; y = scalb(y,-k); x=scalb(x,-k); }
 
     /* begin argument reduction */
 begin:
-	if (t < 2.4375) {		 
+	if (t < 2.4375) {
 
 	/* truncate 4(t+1/16) to integer for branching */
 	    k = 4 * (t+0.0625);
 	    switch (k) {
 
 	    /* t is in [0,7/16] */
-	    case 0:                    
+	    case 0:
 	    case 1:
-		if (t < small) 
+		if (t < small)
 		    { big + small ;  /* raise inexact flag */
 		      return (copysign((signx>zero)?t:PI-t,signy)); }
 
 		hi = zero;  lo = zero;  break;
 
 	    /* t is in [7/16,11/16] */
-	    case 2:                    
+	    case 2:
 		hi = athfhi; lo = athflo;
 		z = x+x;
 		t = ( (y+y) - x ) / ( z +  y ); break;
 
 	    /* t is in [11/16,19/16] */
-	    case 3:                    
+	    case 3:
 	    case 4:
 		hi = PIo4; lo = zero;
 		t = ( y - x ) / ( x + y ); break;
 
 	    /* t is in [19/16,39/16] */
-	    default:                   
+	    default:
 		hi = at1fhi; lo = at1flo;
 		z = y-x; y=y+y+y; t = x+x;
 		t = ( (z+z)-x ) / ( t + y ); break;
@@ -252,7 +252,7 @@ begin:
 	}
 	/* end of if (t < 2.4375) */
 
-	else                           
+	else
 	{
 	    hi = PIo2; lo = zero;
 
@@ -260,7 +260,7 @@ begin:
 	    if (t <= big)  t = - x / y;
 
 	    /* t is in [big, INF] */
-	    else          
+	    else
 	      { big+small;	/* raise inexact flag */
 		t = zero; }
 	}
diff --git a/lib/libm/common/sincos.c b/lib/libm/common/sincos.c
index ab88560..fc35618 100644
--- a/lib/libm/common/sincos.c
+++ b/lib/libm/common/sincos.c
@@ -67,7 +67,7 @@ double x;
 }
 
 double
-cos(x) 
+cos(x)
 double x;
 {
 	double a,c,z,s = 1.0;
@@ -83,7 +83,7 @@ double x;
 		}
 		else {			/* ... in [PI/4,3PI/4] */
 			a = PIo2-a;
-			return a+a*sin__S(a*a);	/* rtn. S(PI/2-|x|) */ 
+			return a+a*sin__S(a*a);	/* rtn. S(PI/2-|x|) */
 		}
 	}
 	if (a < small) {
diff --git a/lib/libm/common/tan.c b/lib/libm/common/tan.c
index 61ed5c5..7b49bce 100644
--- a/lib/libm/common/tan.c
+++ b/lib/libm/common/tan.c
@@ -37,7 +37,7 @@ static char sccsid[] = "@(#)tan.c	8.1 (Berkeley) 6/4/93";
 
 #include "trig.h"
 double
-tan(x) 
+tan(x)
 double x;
 {
 	double a,z,ss,cc,c;
diff --git a/lib/libm/common/trig.h b/lib/libm/common/trig.h
index 9e05b0e..e31fb4c 100644
--- a/lib/libm/common/trig.h
+++ b/lib/libm/common/trig.h
@@ -67,7 +67,7 @@ static const double
 	zero = 0,
 	one = 1,
 	negone = -1,
-	half = 1.0/2.0, 
+	half = 1.0/2.0,
 	small = 1E-10,	/* 1+small**2 == 1; better values for small:
 			 *		small	= 1.5E-9 for VAX D
 			 *			= 1.2E-8 for IEEE Double
@@ -77,27 +77,27 @@ static const double
 
 /* sin__S(x*x) ... re-implemented as a macro
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) 
- * CODED IN C BY K.C. NG, 1/21/85; 
+ * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X)
+ * CODED IN C BY K.C. NG, 1/21/85;
  * REVISED BY K.C. NG on 8/13/85.
  *
  *	    sin(x*k) - x
  * RETURN  --------------- on [-PI/4,PI/4] , where k=pi/PI, PI is the rounded
- *	            x	
+ *	            x
  * value of pi in machine precision:
  *
  *	Decimal:
- *		pi = 3.141592653589793 23846264338327 ..... 
+ *		pi = 3.141592653589793 23846264338327 .....
  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
+ *    56 bits   PI = 3.141592653589793 227020265 ..... ,
  *
  *	Hexadecimal:
  *		pi = 3.243F6A8885A308D313198A2E....
  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18
- *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    
+ *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2
  *
  * Method:
- *	1. Let z=x*x. Create a polynomial approximation to 
+ *	1. Let z=x*x. Create a polynomial approximation to
  *	    (sin(k*x)-x)/x  =  z*(S0 + S1*z^1 + ... + S5*z^5).
  *	Then
  *      sin__S(x*x) = z*(S0 + S1*z^1 + ... + S5*z^5)
@@ -105,8 +105,8 @@ static const double
  *	The coefficient S's are obtained by a special Remez algorithm.
  *
  * Accuracy:
- *	In the absence of rounding error, the approximation has absolute error 
- *	less than 2**(-61.11) for VAX D FORMAT, 2**(-57.45) for IEEE DOUBLE. 
+ *	In the absence of rounding error, the approximation has absolute error
+ *	less than 2**(-61.11) for VAX D FORMAT, 2**(-57.45) for IEEE DOUBLE.
  *
  * Constants:
  * The hexadecimal values are the intended ones for the following constants.
@@ -149,28 +149,28 @@ ic(S5,  1.5868926979889205164E-10 , -33,  1.5CF61DF672B13)
 
 /* cos__C(x*x) ... re-implemented as a macro
  * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
- * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X) 
- * CODED IN C BY K.C. NG, 1/21/85; 
+ * STATIC KERNEL FUNCTION OF SIN(X), COS(X), AND TAN(X)
+ * CODED IN C BY K.C. NG, 1/21/85;
  * REVISED BY K.C. NG on 8/13/85.
  *
- *	   		    x*x	
+ *	   		    x*x
  * RETURN   cos(k*x) - 1 + ----- on [-PI/4,PI/4],  where k = pi/PI,
- *	  		     2	
+ *	  		     2
  * PI is the rounded value of pi in machine precision :
  *
  *	Decimal:
- *		pi = 3.141592653589793 23846264338327 ..... 
+ *		pi = 3.141592653589793 23846264338327 .....
  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
+ *    56 bits   PI = 3.141592653589793 227020265 ..... ,
  *
  *	Hexadecimal:
  *		pi = 3.243F6A8885A308D313198A2E....
  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18
- *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    
+ *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2
  *
  *
  * Method:
- *	1. Let z=x*x. Create a polynomial approximation to 
+ *	1. Let z=x*x. Create a polynomial approximation to
  *	    cos(k*x)-1+z/2  =  z*z*(C0 + C1*z^1 + ... + C5*z^5)
  *	then
  *      cos__C(z) =  z*z*(C0 + C1*z^1 + ... + C5*z^5)
@@ -178,9 +178,9 @@ ic(S5,  1.5868926979889205164E-10 , -33,  1.5CF61DF672B13)
  *	The coefficient C's are obtained by a special Remez algorithm.
  *
  * Accuracy:
- *	In the absence of rounding error, the approximation has absolute error 
- *	less than 2**(-64) for VAX D FORMAT, 2**(-58.3) for IEEE DOUBLE. 
- *	
+ *	In the absence of rounding error, the approximation has absolute error
+ *	less than 2**(-64) for VAX D FORMAT, 2**(-58.3) for IEEE DOUBLE.
+ *
  *
  * Constants:
  * The hexadecimal values are the intended ones for the following constants.
diff --git a/lib/libm/common_source/acosh.c b/lib/libm/common_source/acosh.c
index bc16cc7..149e5de 100644
--- a/lib/libm/common_source/acosh.c
+++ b/lib/libm/common_source/acosh.c
@@ -48,10 +48,10 @@ static char sccsid[] = "@(#)acosh.c	8.1 (Berkeley) 6/4/93";
  *	log1p(x) 		...return log(1+x)
  *
  * Method :
- *	Based on 
+ *	Based on
  *		acosh(x) = log [ x + sqrt(x*x-1) ]
  *	we have
- *		acosh(x) := log1p(x)+ln2,	if (x > 1.0E20); else		
+ *		acosh(x) := log1p(x)+ln2,	if (x > 1.0E20); else
  *		acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) .
  *	These formulae avoid the over/underflow complication.
  *
@@ -60,7 +60,7 @@ static char sccsid[] = "@(#)acosh.c	8.1 (Berkeley) 6/4/93";
  *	acosh(NaN) is NaN without signal.
  *
  * Accuracy:
- *	acosh(x) returns the exact inverse hyperbolic cosine of x nearly 
+ *	acosh(x) returns the exact inverse hyperbolic cosine of x nearly
  *	rounded. In a test run with 512,000 random arguments on a VAX, the
  *	maximum observed error was 3.30 ulps (units of the last place) at
  *	x=1.0070493753568216 .
@@ -87,7 +87,7 @@ ic(ln2lo, 1.9082149292705877000E-10,-33, 1.A39EF35793C76)
 
 double acosh(x)
 double x;
-{	
+{
 	double t,big=1.E20; /* big+1==big */
 
 #if !defined(vax)&&!defined(tahoe)
@@ -95,7 +95,7 @@ double x;
 #endif	/* !defined(vax)&&!defined(tahoe) */
 
     /* return log1p(x) + log(2) if x is large */
-	if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);} 
+	if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);}
 
 	t=sqrt(x-1.0);
 	return(log1p(t*(t+sqrt(x+1.0))));
diff --git a/lib/libm/common_source/asincos.c b/lib/libm/common_source/asincos.c
index c746b16..12d0e14 100644
--- a/lib/libm/common_source/asincos.c
+++ b/lib/libm/common_source/asincos.c
@@ -45,12 +45,12 @@ static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";
  *	sqrt(x)
  *
  * Required kernel function:
- *	atan2(y,x) 
+ *	atan2(y,x)
  *
- * Method :                  
- *	asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is 
+ * Method :
+ *	asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is
  *		  computed as follows
- *			1-x*x                     if x <  0.5, 
+ *			1-x*x                     if x <  0.5,
  *			2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.
  *
  * Special cases:
@@ -59,22 +59,22 @@ static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";
  *
  * Accuracy:
  * 1)  If atan2() uses machine PI, then
- * 
+ *
  *	asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;
  *	and PI is the exact pi rounded to machine precision (see atan2 for
  *      details):
  *
  *	in decimal:
- *		pi = 3.141592653589793 23846264338327 ..... 
+ *		pi = 3.141592653589793 23846264338327 .....
  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
+ *    56 bits   PI = 3.141592653589793 227020265 ..... ,
  *
  *	in hexadecimal:
  *		pi = 3.243F6A8885A308D313198A2E....
  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *	
- *	In a test run with more than 200,000 random arguments on a VAX, the 
+ *
+ *	In a test run with more than 200,000 random arguments on a VAX, the
  *	maximum observed error in ulps (units in the last place) was
  *	2.06 ulps.      (comparing against (PI/pi)*(exact asin(x)));
  *
@@ -82,7 +82,7 @@ static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";
  *
  *	asin(x) returns the exact asin(x) with error below about 2 ulps.
  *
- *	In a test run with more than 1,024,000 random arguments on a VAX, the 
+ *	In a test run with more than 1,024,000 random arguments on a VAX, the
  *	maximum observed error in ulps (units in the last place) was
  *      1.99 ulps.
  */
@@ -97,7 +97,7 @@ double x;
 	s=copysign(x,one);
 	if(s <= 0.5)
 	    return(atan2(x,sqrt(one-x*x)));
-	else 
+	else
 	    { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }
 
 }
@@ -112,9 +112,9 @@ double x;
  *	sqrt(x)
  *
  * Required kernel function:
- *	atan2(y,x) 
+ *	atan2(y,x)
  *
- * Method :                  
+ * Method :
  *			      ________
  *                           / 1 - x
  *	acos(x) = 2*atan2(  / -------- , 1 ) .
@@ -126,22 +126,22 @@ double x;
  *
  * Accuracy:
  * 1)  If atan2() uses machine PI, then
- * 
+ *
  *	acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;
  *	and PI is the exact pi rounded to machine precision (see atan2 for
  *      details):
  *
  *	in decimal:
- *		pi = 3.141592653589793 23846264338327 ..... 
+ *		pi = 3.141592653589793 23846264338327 .....
  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
+ *    56 bits   PI = 3.141592653589793 227020265 ..... ,
  *
  *	in hexadecimal:
  *		pi = 3.243F6A8885A308D313198A2E....
  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *	
- *	In a test run with more than 200,000 random arguments on a VAX, the 
+ *
+ *	In a test run with more than 200,000 random arguments on a VAX, the
  *	maximum observed error in ulps (units in the last place) was
  *	2.07 ulps.      (comparing against (PI/pi)*(exact acos(x)));
  *
@@ -149,7 +149,7 @@ double x;
  *
  *	acos(x) returns the exact acos(x) with error below about 2 ulps.
  *
- *	In a test run with more than 1,024,000 random arguments on a VAX, the 
+ *	In a test run with more than 1,024,000 random arguments on a VAX, the
  *	maximum observed error in ulps (units in the last place) was
  *	2.15 ulps.
  */
diff --git a/lib/libm/common_source/asinh.c b/lib/libm/common_source/asinh.c
index 5db8d2d..1804145 100644
--- a/lib/libm/common_source/asinh.c
+++ b/lib/libm/common_source/asinh.c
@@ -49,16 +49,16 @@ static char sccsid[] = "@(#)asinh.c	8.1 (Berkeley) 6/4/93";
  *	log1p(x) 		...return log(1+x)
  *
  * Method :
- *	Based on 
+ *	Based on
  *		asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
  *	we have
  *	asinh(x) := x  if  1+x*x=1,
  *		 := sign(x)*(log1p(x)+ln2))	 if sqrt(1+x*x)=x, else
- *		 := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) )  
+ *		 := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) )
  *
  * Accuracy:
  *	asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded.
- *	In a test run with 52,000 random arguments on a VAX, the maximum 
+ *	In a test run with 52,000 random arguments on a VAX, the maximum
  *	observed error was 1.58 ulps (units in the last place).
  *
  * Constants:
@@ -82,16 +82,16 @@ ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76)
 
 double asinh(x)
 double x;
-{	
+{
 	double t,s;
 	const static double	small=1.0E-10,	/* fl(1+small*small) == 1 */
 				big  =1.0E20,	/* fl(1+big) == big */
-				one  =1.0   ;	
+				one  =1.0   ;
 
 #if !defined(vax)&&!defined(tahoe)
 	if(x!=x) return(x);	/* x is NaN */
 #endif	/* !defined(vax)&&!defined(tahoe) */
-	if((t=copysign(x,one))>small) 
+	if((t=copysign(x,one))>small)
 	    if(t<big) {
 	     	s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); }
 	    else	/* if |x| > big */
diff --git a/lib/libm/common_source/atan.c b/lib/libm/common_source/atan.c
index 272c7f1..f29c7d4 100644
--- a/lib/libm/common_source/atan.c
+++ b/lib/libm/common_source/atan.c
@@ -41,32 +41,32 @@ static char sccsid[] = "@(#)atan.c	8.1 (Berkeley) 6/4/93";
  * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.
  *
  * Required kernel function:
- *	atan2(y,x) 
+ *	atan2(y,x)
  *
- * Method:                  
- *	atan(x) = atan2(x,1.0). 
+ * Method:
+ *	atan(x) = atan2(x,1.0).
  *
  * Special case:
  *	if x is NaN, return x itself.
  *
  * Accuracy:
  * 1)  If atan2() uses machine PI, then
- * 
+ *
  *	atan(x) returns (PI/pi) * (the exact arc tangent of x) nearly rounded;
  *	and PI is the exact pi rounded to machine precision (see atan2 for
  *      details):
  *
  *	in decimal:
- *		pi = 3.141592653589793 23846264338327 ..... 
+ *		pi = 3.141592653589793 23846264338327 .....
  *    53 bits   PI = 3.141592653589793 115997963 ..... ,
- *    56 bits   PI = 3.141592653589793 227020265 ..... ,  
+ *    56 bits   PI = 3.141592653589793 227020265 ..... ,
  *
  *	in hexadecimal:
  *		pi = 3.243F6A8885A308D313198A2E....
  *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps
  *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps
- *	
- *	In a test run with more than 200,000 random arguments on a VAX, the 
+ *
+ *	In a test run with more than 200,000 random arguments on a VAX, the
  *	maximum observed error in ulps (units in the last place) was
  *	0.86 ulps.      (comparing against (PI/pi)*(exact atan(x))).
  *
@@ -74,7 +74,7 @@ static char sccsid[] = "@(#)atan.c	8.1 (Berkeley) 6/4/93";
  *
  *	atan(x) returns the exact atan(x) with error below about 2 ulps.
  *
- *	In a test run with more than 1,024,000 random arguments on a VAX, the 
+ *	In a test run with more than 1,024,000 random arguments on a VAX, the
  *	maximum observed error in ulps (units in the last place) was
  *	0.85 ulps.
  */
diff --git a/lib/libm/common_source/atanh.c b/lib/libm/common_source/atanh.c
index 89cb61c..e5cdadd 100644
--- a/lib/libm/common_source/atanh.c
+++ b/lib/libm/common_source/atanh.c
@@ -38,14 +38,14 @@ static char sccsid[] = "@(#)atanh.c	8.1 (Berkeley) 6/4/93";
 /* ATANH(X)
  * RETURN THE HYPERBOLIC ARC TANGENT OF X
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85; 
+ * CODED IN C BY K.C. NG, 1/8/85;
  * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85.
  *
  * Required kernel function:
  *	log1p(x) 	...return log(1+x)
  *
  * Method :
- *	Return 
+ *	Return
  *                          1              2x                          x
  *		atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  *                          2             1 - x                      1 - x
diff --git a/lib/libm/common_source/cosh.c b/lib/libm/common_source/cosh.c
index e2b3073..e8d3519 100644
--- a/lib/libm/common_source/cosh.c
+++ b/lib/libm/common_source/cosh.c
@@ -38,7 +38,7 @@ static char sccsid[] = "@(#)cosh.c	8.1 (Berkeley) 6/4/93";
 /* COSH(X)
  * RETURN THE HYPERBOLIC COSINE OF X
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85; 
+ * CODED IN C BY K.C. NG, 1/8/85;
  * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85.
  *
  * Required system supported functions :
@@ -46,20 +46,20 @@ static char sccsid[] = "@(#)cosh.c	8.1 (Berkeley) 6/4/93";
  *	scalb(x,N)
  *
  * Required kernel function:
- *	exp(x) 
+ *	exp(x)
  *	exp__E(x,c)	...return exp(x+c)-1-x for |x|<0.3465
  *
  * Method :
- *	1. Replace x by |x|. 
- *	2. 
- *		                                        [ exp(x) - 1 ]^2 
+ *	1. Replace x by |x|.
+ *	2.
+ *		                                        [ exp(x) - 1 ]^2
  *	    0        <= x <= 0.3465  :  cosh(x) := 1 + -------------------
  *			       			           2*exp(x)
  *
  *		                                   exp(x) +  1/exp(x)
  *	    0.3465   <= x <= 22      :  cosh(x) := -------------------
  *			       			           2
- *	    22       <= x <= lnovfl  :  cosh(x) := exp(x)/2 
+ *	    22       <= x <= lnovfl  :  cosh(x) := exp(x)/2
  *	    lnovfl   <= x <= lnovfl+log(2)
  *				     :  cosh(x) := exp(x)/2 (avoid overflow)
  *	    log(2)+lnovfl <  x <  INF:  overflow to INF
@@ -106,7 +106,7 @@ static max = 1023                     ;
 
 double cosh(x)
 double x;
-{	
+{
 	static const double half=1.0/2.0,
 		one=1.0, small=1.0E-18; /* fl(1+small)==1 */
 	double t;
@@ -115,19 +115,19 @@ double x;
 	if(x!=x) return(x);	/* x is NaN */
 #endif	/* !defined(vax)&&!defined(tahoe) */
 	if((x=copysign(x,one)) <= 22)
-	    if(x<0.3465) 
+	    if(x<0.3465)
 		if(x<small) return(one+x);
 		else {t=x+__exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); }
 
 	    else /* for x lies in [0.3465,22] */
 	        { t=exp(x); return((t+one/t)*half); }
 
-	if( lnovfl <= x && x <= (lnovfl+0.7)) 
-        /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1)) 
-         * and return 2^max*exp(x) to avoid unnecessary overflow 
+	if( lnovfl <= x && x <= (lnovfl+0.7))
+        /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1))
+         * and return 2^max*exp(x) to avoid unnecessary overflow
          */
-	    return(scalb(exp((x-mln2hi)-mln2lo), max)); 
+	    return(scalb(exp((x-mln2hi)-mln2lo), max));
 
-	else 
+	else
 	    return(exp(x)*half);	/* for large x,  cosh(x)=exp(x)/2 */
 }
diff --git a/lib/libm/common_source/erf.c b/lib/libm/common_source/erf.c
index 308f1a9..ba0c006 100644
--- a/lib/libm/common_source/erf.c
+++ b/lib/libm/common_source/erf.c
@@ -68,13 +68,13 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";
  *          erfc(x) = 1 - erf(x)           if x<=0.25
  *                  = 0.5 + ((0.5-x)-x*P)  if x in [0.25,0.84375]
  *	   where
- *			2		 2	  4		  20  
+ *			2		 2	  4		  20
  *              P =  P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x  )
  * 	   is an approximation to (erf(x)-x)/x with precision
  *
  *						 -56.45
  *			| P - (erf(x)-x)/x | <= 2
- *	
+ *
  *
  *	   Remark. The formula is derived by noting
  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
@@ -96,7 +96,7 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";
  *	   That is, we use rational approximation to approximate
  *			erf(1+s) - (c = (single)0.84506291151)
  *	   Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
- *	   where 
+ *	   where
  *		P1(s) = degree 6 poly in s
  *		Q1(s) = degree 6 poly in s
  *
@@ -105,7 +105,7 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";
  *		erfc(x)	= (1/x)exp(-x*x-(.5*log(pi) -.5z + R(z)/S(z))
  *
  *	Where z = 1/(x*x), R is degree 9, and S is degree 3;
- *	
+ *
  *      5. For x in [4,28]
  *         	erf(x)  = 1.0 - tiny
  *		erfc(x)	= (1/x)exp(-x*x-(.5*log(pi)+eps + zP(z))
@@ -139,7 +139,7 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";
  *
  *      7. Special cases:
  *         	erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
- *         	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, 
+ *         	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
  *	   	erfc/erf(NaN) is NaN
  */
 
@@ -181,7 +181,7 @@ p8 =   1.640186161764254363152286358441771740838e-0006,
 p9 =  -1.571599331700515057841960987689515895479e-0007,
 p10=   1.073087585213621540635426191486561494058e-0008;
 /*
- * Coefficients for approximation to erf in [0.84375,1.25] 
+ * Coefficients for approximation to erf in [0.84375,1.25]
  */
 static double
 pa0 =  -2.362118560752659485957248365514511540287e-0003,
@@ -319,7 +319,7 @@ double erf(x)
 		return (z-one);
 }
 
-double erfc(x) 
+double erfc(x)
 	double x;
 {
 	double R,S,P,Q,s,ax,y,z,r,fabs(),__exp__D();
@@ -352,7 +352,7 @@ double erfc(x)
 	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
 	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
 	    if (x>=0) {
-	        z  = one-c; return z - P/Q; 
+	        z  = one-c; return z - P/Q;
 	    } else {
 		z = c+P/Q; return one+z;
 	    }
diff --git a/lib/libm/common_source/exp.c b/lib/libm/common_source/exp.c
index 9b4f045..fa6ea75 100644
--- a/lib/libm/common_source/exp.c
+++ b/lib/libm/common_source/exp.c
@@ -38,21 +38,21 @@ static char sccsid[] = "@(#)exp.c	8.1 (Berkeley) 6/4/93";
 /* EXP(X)
  * RETURN THE EXPONENTIAL OF X
  * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
- * CODED IN C BY K.C. NG, 1/19/85; 
+ * CODED IN C BY K.C. NG, 1/19/85;
  * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
  *
  * Required system supported functions:
- *	scalb(x,n)	
- *	copysign(x,y)	
+ *	scalb(x,n)
+ *	copysign(x,y)
  *	finite(x)
  *
  * Method:
- *	1. Argument Reduction: given the input x, find r and integer k such 
+ *	1. Argument Reduction: given the input x, find r and integer k such
  *	   that
- *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .  
+ *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .
  *	   r will be represented as r := z+c for better accuracy.
  *
- *	2. Compute exp(r) by 
+ *	2. Compute exp(r) by
  *
  *		exp(r) = 1 + r + r*R1/(2-R1),
  *	   where
@@ -143,7 +143,7 @@ double x;
 		}
 		/* end of x > lntiny */
 
-		else 
+		else
 		     /* exp(-big#) underflows to zero */
 		     if(finite(x))  return(scalb(1.0,-5000));
 
@@ -152,7 +152,7 @@ double x;
 	}
 	/* end of x < lnhuge */
 
-	else 
+	else
 	/* exp(INF) is INF, exp(+big#) overflows to INF */
 	    return( finite(x) ?  scalb(1.0,5000)  : x);
 }
@@ -188,7 +188,7 @@ double x, c;
 		}
 		/* end of x > lntiny */
 
-		else 
+		else
 		     /* exp(-big#) underflows to zero */
 		     if(finite(x))  return(scalb(1.0,-5000));
 
@@ -197,7 +197,7 @@ double x, c;
 	}
 	/* end of x < lnhuge */
 
-	else 
+	else
 	/* exp(INF) is INF, exp(+big#) overflows to INF */
 	    return( finite(x) ?  scalb(1.0,5000)  : x);
 }
diff --git a/lib/libm/common_source/exp__E.c b/lib/libm/common_source/exp__E.c
index ab97247..16ec4cbb 100644
--- a/lib/libm/common_source/exp__E.c
+++ b/lib/libm/common_source/exp__E.c
@@ -41,7 +41,7 @@ static char sccsid[] = "@(#)exp__E.c	8.1 (Berkeley) 6/4/93";
  * exp__E RETURNS
  *
  *			 /  exp(x+c) - 1 - x ,  1E-19 < |x| < .3465736
- *       exp__E(x,c) = 	| 		     
+ *       exp__E(x,c) = 	|
  *			 \  0 ,  |x| < 1E-19.
  *
  * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
@@ -50,12 +50,12 @@ static char sccsid[] = "@(#)exp__E.c	8.1 (Berkeley) 6/4/93";
  * REVISED BY K.C. NG on 3/16/85, 4/16/85.
  *
  * Required system supported function:
- *	copysign(x,y)	
+ *	copysign(x,y)
  *
  * Method:
  *	1. Rational approximation. Let r=x+c.
  *	   Based on
- *                                   2 * sinh(r/2)     
+ *                                   2 * sinh(r/2)
  *                exp(r) - 1 =   ----------------------   ,
  *                               cosh(r/2) - sinh(r/2)
  *	   exp__E(r) is computed using
@@ -76,7 +76,7 @@ static char sccsid[] = "@(#)exp__E.c	8.1 (Berkeley) 6/4/93";
  * Approximation error:
  *
  *   |	exp(x) - 1			   |        2**(-57),  (IEEE double)
- *   | ------------  -  (exp__E(x,0)+x)/x  |  <= 
+ *   | ------------  -  (exp__E(x,0)+x)/x  |  <=
  *   |	     x			           |	    2**(-69).  (VAX D)
  *
  * Constants:
@@ -120,7 +120,7 @@ double x,c;
 #else	/* defined(vax)||defined(tahoe) */
            q = z*( q1 +z*  q2 );
 #endif	/* defined(vax)||defined(tahoe) */
-           xp= x*p     ; 
+           xp= x*p     ;
 	   xh= x*half  ;
            w = xh-(q-xp)  ;
 	   p = p+p;
diff --git a/lib/libm/common_source/expm1.c b/lib/libm/common_source/expm1.c
index 760d2be..8a84f14 100644
--- a/lib/libm/common_source/expm1.c
+++ b/lib/libm/common_source/expm1.c
@@ -38,36 +38,36 @@ static char sccsid[] = "@(#)expm1.c	8.1 (Berkeley) 6/4/93";
 /* EXPM1(X)
  * RETURN THE EXPONENTIAL OF X MINUS ONE
  * DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 BITS)
- * CODED IN C BY K.C. NG, 1/19/85; 
+ * CODED IN C BY K.C. NG, 1/19/85;
  * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/21/85, 4/16/85.
  *
  * Required system supported functions:
- *	scalb(x,n)	
- *	copysign(x,y)	
+ *	scalb(x,n)
+ *	copysign(x,y)
  *	finite(x)
  *
  * Kernel function:
  *	exp__E(x,c)
  *
  * Method:
- *	1. Argument Reduction: given the input x, find r and integer k such 
+ *	1. Argument Reduction: given the input x, find r and integer k such
  *	   that
- *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .  
+ *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .
  *	   r will be represented as r := z+c for better accuracy.
  *
- *	2. Compute EXPM1(r)=exp(r)-1 by 
+ *	2. Compute EXPM1(r)=exp(r)-1 by
  *
  *			EXPM1(r=z+c) := z + exp__E(z,c)
  *
  *	3. EXPM1(x) =  2^k * ( EXPM1(r) + 1-2^-k ).
  *
- * 	Remarks: 
+ * 	Remarks:
  *	   1. When k=1 and z < -0.25, we use the following formula for
  *	      better accuracy:
  *			EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) )
  *	   2. To avoid rounding error in 1-2^-k where k is large, we use
  *			EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 }
- *	      when k>56. 
+ *	      when k>56.
  *
  * Special cases:
  *	EXPM1(INF) is INF, EXPM1(NaN) is NaN;
@@ -108,7 +108,7 @@ ic(invln2, 1.4426950408889633870E0,     0, 1.71547652B82FE)
 double expm1(x)
 double x;
 {
-	const static double one=1.0, half=1.0/2.0; 
+	const static double one=1.0, half=1.0/2.0;
 	double  z,hi,lo,c;
 	int k;
 #if defined(vax)||defined(tahoe)
@@ -126,13 +126,13 @@ double x;
 
 		    /* argument reduction : x - k*ln2 */
 			k= invln2 *x+copysign(0.5,x);	/* k=NINT(x/ln2) */
-			hi=x-k*ln2hi ; 
+			hi=x-k*ln2hi ;
 			z=hi-(lo=k*ln2lo);
 			c=(hi-z)-lo;
 
 			if(k==0) return(z+__exp__E(z,c));
 			if(k==1)
-			    if(z< -0.25) 
+			    if(z< -0.25)
 				{x=z+half;x +=__exp__E(z,c); return(x+x);}
 			    else
 				{z+=__exp__E(z,c); x=half+z; return(x+x);}
@@ -143,17 +143,17 @@ double x;
 			      { x=one-scalb(one,-k); z += __exp__E(z,c);}
 			    else if(k<100)
 			      { x = __exp__E(z,c)-scalb(one,-k); x+=z; z=one;}
-			    else 
+			    else
 			      { x = __exp__E(z,c)+z; z=one;}
 
-			    return (scalb(x+z,k));  
+			    return (scalb(x+z,k));
 			}
 		}
 		/* end of x > lnunfl */
 
-		else 
+		else
 		     /* expm1(-big#) rounded to -1 (inexact) */
-		     if(finite(x))  
+		     if(finite(x))
 			 { ln2hi+ln2lo; return(-one);}
 
 		     /* expm1(-INF) is -1 */
@@ -161,7 +161,7 @@ double x;
 	}
 	/* end of x < lnhuge */
 
-	else 
+	else
 	/*  expm1(INF) is INF, expm1(+big#) overflows to INF */
 	    return( finite(x) ?  scalb(one,5000) : x);
 }
diff --git a/lib/libm/common_source/j0.c b/lib/libm/common_source/j0.c
index a9b28b3..684fb4d 100644
--- a/lib/libm/common_source/j0.c
+++ b/lib/libm/common_source/j0.c
@@ -46,18 +46,18 @@ static char sccsid[] = "@(#)j0.c	8.2 (Berkeley) 11/30/93";
  *
  * 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.
  * ====================================================
  *
  * ******************* WARNING ********************
  * This is an alpha version of SunPro's FDLIBM (Freely
- * Distributable Math Library) for IEEE double precision 
+ * Distributable Math Library) for IEEE double precision
  * arithmetic. FDLIBM is a basic math library written
- * in C that runs on machines that conform to IEEE 
- * Standard 754/854. This alpha version is distributed 
- * for testing purpose. Those who use this software 
- * should report any bugs to 
+ * in C that runs on machines that conform to IEEE
+ * Standard 754/854. This alpha version is distributed
+ * for testing purpose. Those who use this software
+ * should report any bugs to
  *
  *		fdlibm-comments@sunpro.eng.sun.com
  *
@@ -84,20 +84,20 @@ static char sccsid[] = "@(#)j0.c	8.2 (Berkeley) 11/30/93";
  * 	   (To avoid cancellation, use
  *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  * 	    to compute the worse one.)
- *	   
+ *
  *	3 Special cases
  *		j0(nan)= nan
  *		j0(0) = 1
  *		j0(inf) = 0
- *		
+ *
  * Method -- y0(x):
  *	1. For x<2.
- *	   Since 
+ *	   Since
  *		y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
  *	   therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
  *	   We use the following function to approximate y0,
  *		y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
- *	   where 
+ *	   where
  *		U(z) = u0 + u1*z + ... + u6*z^6
  *		V(z) = 1  + v1*z + ... + v4*z^4
  *	   with absolute approximation error bounded by 2**-72.
@@ -121,7 +121,7 @@ static char sccsid[] = "@(#)j0.c	8.2 (Berkeley) 11/30/93";
 
 static double pzero __P((double)), qzero __P((double));
 
-static double 
+static double
 huge 	= 1e300,
 zero    = 0.0,
 one	= 1.0,
@@ -138,7 +138,7 @@ s03 =   5.135465502073181376284426245689510134134e-0007,
 s04 =   1.166140033337900097836930825478674320464e-0009;
 
 double
-j0(x) 
+j0(x)
 	double x;
 {
 	double z, s,c,ss,cc,r,u,v;
@@ -201,7 +201,7 @@ v03 =   2.591508518404578033173189144579208685163e-0007,
 v04 =   4.411103113326754838596529339004302243157e-0010;
 
 double
-y0(x) 
+y0(x)
 	double x;
 {
 	double z, s, c, ss, cc, u, v;
@@ -345,7 +345,7 @@ static double pzero(x)
 	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.
diff --git a/lib/libm/common_source/j1.c b/lib/libm/common_source/j1.c
index 71602aa..e8ca43a 100644
--- a/lib/libm/common_source/j1.c
+++ b/lib/libm/common_source/j1.c
@@ -46,18 +46,18 @@ static char sccsid[] = "@(#)j1.c	8.2 (Berkeley) 11/30/93";
  *
  * 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.
  * ====================================================
  *
  * ******************* WARNING ********************
  * This is an alpha version of SunPro's FDLIBM (Freely
- * Distributable Math Library) for IEEE double precision 
+ * Distributable Math Library) for IEEE double precision
  * arithmetic. FDLIBM is a basic math library written
- * in C that runs on machines that conform to IEEE 
- * Standard 754/854. This alpha version is distributed 
- * for testing purpose. Those who use this software 
- * should report any bugs to 
+ * in C that runs on machines that conform to IEEE
+ * Standard 754/854. This alpha version is distributed
+ * for testing purpose. Those who use this software
+ * should report any bugs to
  *
  *		fdlibm-comments@sunpro.eng.sun.com
  *
@@ -85,16 +85,16 @@ static char sccsid[] = "@(#)j1.c	8.2 (Berkeley) 11/30/93";
  * 	   (To avoid cancellation, use
  *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  * 	    to compute the worse one.)
- *	   
+ *
  *	3 Special cases
  *		j1(nan)= nan
  *		j1(0) = 0
  *		j1(inf) = 0
- *		
+ *
  * Method -- y1(x):
- *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN 
+ *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
  *	2. For x<2.
- *	   Since 
+ *	   Since
  *		y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
  *	   therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
  *	   We use the following function to approximate y1,
@@ -122,7 +122,7 @@ static char sccsid[] = "@(#)j1.c	8.2 (Berkeley) 11/30/93";
 
 static double pone(), qone();
 
-static double 
+static double
 huge    = 1e300,
 zero    = 0.0,
 one	= 1.0,
@@ -142,7 +142,7 @@ s05 =   1.235422744261379203512624973117299248281e-0011;
 
 #define two_129	6.80564733841876926e+038	/* 2^129 */
 #define two_m54	5.55111512312578270e-017	/* 2^-54 */
-double j1(x) 
+double j1(x)
 	double x;
 {
 	double z, s,c,ss,cc,r,u,v,y;
@@ -205,7 +205,7 @@ static double v0[5] = {
    1.665592462079920695971450872592458916421e-0011,
 };
 
-double y1(x) 
+double y1(x)
 	double x;
 {
 	double z, s, c, ss, cc, u, v;
@@ -254,10 +254,10 @@ double y1(x)
                     z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
                 }
                 return z;
-        } 
+        }
         if (x <= two_m54) {    /* x < 2**-54 */
             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]))));
@@ -351,7 +351,7 @@ static double pone(x)
 	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.
diff --git a/lib/libm/common_source/jn.c b/lib/libm/common_source/jn.c
index 85a5401..28d9687 100644
--- a/lib/libm/common_source/jn.c
+++ b/lib/libm/common_source/jn.c
@@ -46,18 +46,18 @@ static char sccsid[] = "@(#)jn.c	8.2 (Berkeley) 11/30/93";
  *
  * 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.
  * ====================================================
  *
  * ******************* WARNING ********************
  * This is an alpha version of SunPro's FDLIBM (Freely
- * Distributable Math Library) for IEEE double precision 
+ * Distributable Math Library) for IEEE double precision
  * arithmetic. FDLIBM is a basic math library written
- * in C that runs on machines that conform to IEEE 
- * Standard 754/854. This alpha version is distributed 
- * for testing purpose. Those who use this software 
- * should report any bugs to 
+ * in C that runs on machines that conform to IEEE
+ * Standard 754/854. This alpha version is distributed
+ * for testing purpose. Those who use this software
+ * should report any bugs to
  *
  *		fdlibm-comments@sunpro.eng.sun.com
  *
@@ -69,7 +69,7 @@ static char sccsid[] = "@(#)jn.c	8.2 (Berkeley) 11/30/93";
  * jn(int n, double x), yn(int n, double 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.
@@ -88,7 +88,7 @@ static char sccsid[] = "@(#)jn.c	8.2 (Berkeley) 11/30/93";
  *	yn(n,x) is similar in all respects, except
  *	that forward recursion is used for all
  *	values of n>1.
- *	
+ *
  */
 
 #include <math.h>
@@ -120,7 +120,7 @@ double jn(n,x)
      */
     /* if J(n,NaN) is NaN */
 	if (_IEEE && isnan(x)) return x+x;
-	if (n<0){		
+	if (n<0){
 		n = -n;
 		x = -x;
 	}
@@ -134,10 +134,10 @@ double jn(n,x)
 			/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 	    if (_IEEE && x >= 8.148143905337944345e+090) {
 					/* 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
@@ -154,7 +154,7 @@ double jn(n,x)
 		    case 3: temp =  cos(x)-sin(x); break;
 		}
 		b = invsqrtpi*temp/sqrt(x);
-	    } else {	
+	    } else {
 	        a = j0(x);
 	        b = j1(x);
 	        for(i=1;i<n;i++){
@@ -165,7 +165,7 @@ double jn(n,x)
 	    }
 	} else {
 	    if (x < 1.86264514923095703125e-009) { /* 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 */
@@ -180,14 +180,14 @@ double jn(n,x)
 		}
 	    } 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
@@ -203,9 +203,9 @@ double jn(n,x)
 		 * 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;
@@ -254,7 +254,7 @@ double jn(n,x)
 	}
 	return ((sgn == 1) ? -b : b);
 }
-double yn(n,x) 
+double yn(n,x)
 	int n; double x;
 {
 	int i, sign;
@@ -275,10 +275,10 @@ double yn(n,x)
 	if (n == 0) return(y0(x));
 	if (n == 1) return(sign*y1(x));
 	if(_IEEE && x >= 8.148143905337944345e+090) { /* 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
diff --git a/lib/libm/common_source/log.c b/lib/libm/common_source/log.c
index ae18672..908b8544 100644
--- a/lib/libm/common_source/log.c
+++ b/lib/libm/common_source/log.c
@@ -391,7 +391,7 @@ log(x) double x;
 			return (x+x);
 		else
 			return (infnan(ERANGE));
-	
+
 	/* Argument reduction: 1 <= g < 2; x/2^m = g;	*/
 	/* y = F*(1 + f/F) for |f| <= 2^-8		*/
 
diff --git a/lib/libm/common_source/log10.c b/lib/libm/common_source/log10.c
index d2617cc..75205cd 100644
--- a/lib/libm/common_source/log10.c
+++ b/lib/libm/common_source/log10.c
@@ -38,9 +38,9 @@ static char sccsid[] = "@(#)log10.c	8.1 (Berkeley) 6/4/93";
 /* LOG10(X)
  * RETURN THE BASE 10 LOGARITHM OF x
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/20/85; 
+ * CODED IN C BY K.C. NG, 1/20/85;
  * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85.
- * 
+ *
  * Required kernel function:
  *	log(x)
  *
@@ -52,12 +52,12 @@ static char sccsid[] = "@(#)log10.c	8.1 (Berkeley) 6/4/93";
  *    Note:
  *	  [log(10)]   rounded to 56 bits has error  .0895  ulps,
  *	  [1/log(10)] rounded to 53 bits has error  .198   ulps;
- *	  therefore, for better accuracy, in VAX D format, we divide 
- *	  log(x) by log(10), but in IEEE Double format, we multiply 
+ *	  therefore, for better accuracy, in VAX D format, we divide
+ *	  log(x) by log(10), but in IEEE Double format, we multiply
  *	  log(x) by [1/log(10)].
  *
  * Special cases:
- *	log10(x) is NaN with signal if x < 0; 
+ *	log10(x) is NaN with signal if x < 0;
  *	log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
  *	log10(NaN) is that NaN with no signal.
  *
diff --git a/lib/libm/common_source/log1p.c b/lib/libm/common_source/log1p.c
index cbf9fcd..0202667 100644
--- a/lib/libm/common_source/log1p.c
+++ b/lib/libm/common_source/log1p.c
@@ -35,24 +35,24 @@
 static char sccsid[] = "@(#)log1p.c	8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
-/* LOG1P(x) 
+/* LOG1P(x)
  * RETURN THE LOGARITHM OF 1+x
  * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/19/85; 
+ * CODED IN C BY K.C. NG, 1/19/85;
  * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85.
- * 
+ *
  * Required system supported functions:
- *	scalb(x,n) 
+ *	scalb(x,n)
  *	copysign(x,y)
- *	logb(x)	
+ *	logb(x)
  *	finite(x)
  *
  * Required kernel function:
  *	log__L(z)
  *
  * Method :
- *	1. Argument Reduction: find k and f such that 
- *			1+x  = 2^k * (1+f), 
+ *	1. Argument Reduction: find k and f such that
+ *			1+x  = 2^k * (1+f),
  *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
  *
  *	2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
@@ -65,11 +65,11 @@ static char sccsid[] = "@(#)log1p.c	8.1 (Berkeley) 6/4/93";
  *
  *	   See log__L() for the values of the coefficients.
  *
- *	3. Finally,  log(1+x) = k*ln2 + log(1+f).  
+ *	3. Finally,  log(1+x) = k*ln2 + log(1+f).
  *
  *	Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers
- *		   n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last 
- *		   20 bits (for VAX D format), or the last 21 bits ( for IEEE 
+ *		   n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last
+ *		   20 bits (for VAX D format), or the last 21 bits ( for IEEE
  *		   double) is 0. This ensures n*ln2hi is exactly representable.
  *		2. In step 1, f may not be representable. A correction term c
  *	 	   for f is computed. It follows that the correction term for
@@ -83,7 +83,7 @@ static char sccsid[] = "@(#)log1p.c	8.1 (Berkeley) 6/4/93";
  *	only log1p(0)=0 is exact for finite argument.
  *
  * Accuracy:
- *	log1p(x) returns the exact log(1+x) nearly rounded. In a test run 
+ *	log1p(x) returns the exact log(1+x) nearly rounded. In a test run
  *	with 1,536,000 random arguments on a VAX, the maximum observed
  *	error was .846 ulps (units in the last place).
  *
@@ -114,7 +114,7 @@ ic(sqrt2, 1.4142135623730951455E0,     0, 1.6A09E667F3BCD)
 double log1p(x)
 double x;
 {
-	const static double zero=0.0, negone= -1.0, one=1.0, 
+	const static double zero=0.0, negone= -1.0, one=1.0,
 		      half=1.0/2.0, small=1.0E-20;   /* 1+small == 1 */
 	double z,s,t,c;
 	int k;
@@ -129,7 +129,7 @@ double x;
 	   /* argument reduction */
 	      if(copysign(x,one)<small) return(x);
 	      k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k);
-	      if(z+t >= sqrt2 ) 
+	      if(z+t >= sqrt2 )
 		  { k += 1 ; z *= half; t *= half; }
 	      t += negone; x = z + t;
 	      c = (t-x)+z ;		/* correction term for x */
@@ -162,9 +162,9 @@ double x;
     /* end of if (finite(x)) */
 
     /* log(-INF) is NaN */
-	else if(x<0) 
+	else if(x<0)
 	     return(zero/zero);
 
     /* log(+INF) is INF */
-	else return(x);      
+	else return(x);
 }
diff --git a/lib/libm/common_source/log__L.c b/lib/libm/common_source/log__L.c
index c00158f..207cb0d 100644
--- a/lib/libm/common_source/log__L.c
+++ b/lib/libm/common_source/log__L.c
@@ -39,14 +39,14 @@ static char sccsid[] = "@(#)log__L.c	8.1 (Berkeley) 6/4/93";
  *		LOG(1+X) - 2S			       X
  * RETURN      ---------------  WHERE Z = S*S,  S = ------- , 0 <= Z <= .0294...
  *		      S				     2 + X
- *		     
+ *
  * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
  * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS
- * CODED IN C BY K.C. NG, 1/19/85; 
+ * CODED IN C BY K.C. NG, 1/19/85;
  * REVISED BY K.C. Ng, 2/3/85, 4/16/85.
  *
  * Method :
- *	1. Polynomial approximation: let s = x/(2+x). 
+ *	1. Polynomial approximation: let s = x/(2+x).
  *	   Based on log(1+x) = log(1+s) - log(1-s)
  *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  *
@@ -54,11 +54,11 @@ static char sccsid[] = "@(#)log__L.c	8.1 (Berkeley) 6/4/93";
  *
  *	       z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))
  *
- *	   where z=s*s. (See the listing below for Lk's values.) The 
- *	   coefficients are obtained by a special Remez algorithm. 
+ *	   where z=s*s. (See the listing below for Lk's values.) The
+ *	   coefficients are obtained by a special Remez algorithm.
  *
  * Accuracy:
- *	Assuming no rounding error, the maximum magnitude of the approximation 
+ *	Assuming no rounding error, the maximum magnitude of the approximation
  *	error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)
  *	for VAX D format.
  *
diff --git a/lib/libm/common_source/pow.c b/lib/libm/common_source/pow.c
index cb4fc50..9e92985 100644
--- a/lib/libm/common_source/pow.c
+++ b/lib/libm/common_source/pow.c
@@ -35,17 +35,17 @@
 static char sccsid[] = "@(#)pow.c	8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
-/* POW(X,Y)  
- * RETURN X**Y 
+/* POW(X,Y)
+ * RETURN X**Y
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85; 
+ * CODED IN C BY K.C. NG, 1/8/85;
  * REVISED BY K.C. NG on 7/10/85.
  * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92.
  * Required system supported functions:
- *      scalb(x,n)      
- *      logb(x)         
- *	copysign(x,y)	
- *	finite(x)	
+ *      scalb(x,n)
+ *      logb(x)
+ *	copysign(x,y)
+ *	finite(x)
  *	drem(x,y)
  *
  * Required kernel functions:
@@ -56,7 +56,7 @@ static char sccsid[] = "@(#)pow.c	8.1 (Berkeley) 6/4/93";
  *	1. Compute and return log(x) in three pieces:
  *		log(x) = n*ln2 + hi + lo,
  *	   where n is an integer.
- *	2. Perform y*log(x) by simulating muti-precision arithmetic and 
+ *	2. Perform y*log(x) by simulating muti-precision arithmetic and
  *	   return the answer in three pieces:
  *		y*log(x) = m*ln2 + hi + lo,
  *	   where m is an integer.
@@ -94,7 +94,7 @@ static char sccsid[] = "@(#)pow.c	8.1 (Berkeley) 6/4/93";
  *			pow(integer,integer)
  *	always returns the correct integer provided it is representable.
  *	In a test run with 100,000 random arguments with 0 < x, y < 20.0
- *	on a VAX, the maximum observed error was 1.79 ulps (units in the 
+ *	on a VAX, the maximum observed error was 1.79 ulps (units in the
  *	last place).
  *
  * Constants :
@@ -123,7 +123,7 @@ const static double zero=0.0, one=1.0, two=2.0, negone= -1.0;
 
 static double pow_P __P((double, double));
 
-double pow(x,y)  	
+double pow(x,y)
 double x,y;
 {
 	double t;
diff --git a/lib/libm/common_source/sinh.c b/lib/libm/common_source/sinh.c
index 0516849..075172c 100644
--- a/lib/libm/common_source/sinh.c
+++ b/lib/libm/common_source/sinh.c
@@ -38,7 +38,7 @@ static char sccsid[] = "@(#)sinh.c	8.1 (Berkeley) 6/4/93";
 /* SINH(X)
  * RETURN THE HYPERBOLIC SINE OF X
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85; 
+ * CODED IN C BY K.C. NG, 1/8/85;
  * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85.
  *
  * Required system supported functions :
@@ -50,14 +50,14 @@ static char sccsid[] = "@(#)sinh.c	8.1 (Berkeley) 6/4/93";
  *
  * Method :
  *	1. reduce x to non-negative by sinh(-x) = - sinh(x).
- *	2. 
+ *	2.
  *
  *	                                      expm1(x) + expm1(x)/(expm1(x)+1)
  *	    0 <= x <= lnovfl     : sinh(x) := --------------------------------
  *			       		                      2
  *     lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow)
  * lnovfl+ln2 <  x <  INF        :  overflow to INF
- *	
+ *
  *
  * Special cases:
  *	sinh(x) is x if x is +INF, -INF, or NaN.
@@ -112,7 +112,7 @@ double x;
 	    {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));}
 
 	else if(x <= lnovfl+0.7)
-		/* subtract x by ln(2^(max+1)) and return 2^max*exp(x) 
+		/* subtract x by ln(2^(max+1)) and return 2^max*exp(x)
 	    		to avoid unnecessary overflow */
 	    return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign));
 
diff --git a/lib/libm/common_source/tanh.c b/lib/libm/common_source/tanh.c
index d4923b3..6813b55 100644
--- a/lib/libm/common_source/tanh.c
+++ b/lib/libm/common_source/tanh.c
@@ -38,7 +38,7 @@ static char sccsid[] = "@(#)tanh.c	8.1 (Berkeley) 6/4/93";
 /* TANH(X)
  * RETURN THE HYPERBOLIC TANGENT OF X
  * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 1/8/85; 
+ * CODED IN C BY K.C. NG, 1/8/85;
  * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85.
  *
  * Required system supported functions :
@@ -85,7 +85,7 @@ double x;
 
 	sign=copysign(one,x);
 	x=copysign(x,one);
-	if(x < 22.0) 
+	if(x < 22.0)
 	    if( x > one )
 		return(copysign(one-two/(expm1(x+x)+two),sign));
 	    else if ( x > small )
diff --git a/lib/libm/ieee/cabs.c b/lib/libm/ieee/cabs.c
index c4a6fb0..5f1fc62 100644
--- a/lib/libm/ieee/cabs.c
+++ b/lib/libm/ieee/cabs.c
@@ -38,7 +38,7 @@ static char sccsid[] = "@(#)cabs.c	8.1 (Berkeley) 6/4/93";
 /* HYPOT(X,Y)
  * RETURN THE SQUARE ROOT OF X^2 + Y^2  WHERE Z=X+iY
  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
- * CODED IN C BY K.C. NG, 11/28/84; 
+ * CODED IN C BY K.C. NG, 11/28/84;
  * REVISED BY K.C. NG, 7/12/85.
  *
  * Required system supported functions :
@@ -52,16 +52,16 @@ static char sccsid[] = "@(#)cabs.c	8.1 (Berkeley) 6/4/93";
  *	   y if y > x (hence x is never smaller than y).
  *	2. Hypot(x,y) is computed by:
  *	   Case I, x/y > 2
- *		
+ *
  *				       y
  *		hypot = x + -----------------------------
  *			 		    2
  *			    sqrt ( 1 + [x/y]  )  +  x/y
  *
- *	   Case II, x/y <= 2 
+ *	   Case II, x/y <= 2
  *				                   y
  *		hypot = x + --------------------------------------------------
- *				          		     2 
+ *				          		     2
  *				     			[x/y]   -  2
  *			   (sqrt(2)+1) + (x-y)/y + -----------------------------
  *			 		    			  2
@@ -107,7 +107,7 @@ double
 hypot(x,y)
 double x, y;
 {
-	static const double zero=0, one=1, 
+	static const double zero=0, one=1,
 		      small=1.0E-18;	/* fl(1+small)==1 */
 	static const ibig=30;	/* fl(1+2**(2*ibig))==1 */
 	double t,r;
@@ -115,15 +115,15 @@ double x, y;
 
 	if(finite(x))
 	    if(finite(y))
-	    {	
+	    {
 		x=copysign(x,one);
 		y=copysign(y,one);
-		if(y > x) 
+		if(y > x)
 		    { t=x; x=y; y=t; }
 		if(x == zero) return(zero);
 		if(y == zero) return(x);
 		exp= logb(x);
-		if(exp-(int)logb(y) > ibig ) 	
+		if(exp-(int)logb(y) > ibig )
 			/* raise inexact flag and return |x| */
 		   { one+small; return(x); }
 
@@ -144,7 +144,7 @@ double x, y;
 
 	    else if(y==y)   	   /* y is +-INF */
 		     return(copysign(y,one));
-	    else 
+	    else
 		     return(y);	   /* y is NaN and x is finite */
 
 	else if(x==x) 		   /* x is +-INF */
@@ -199,16 +199,16 @@ double x, y;
 
 	if(finite(x))
 	    if(finite(y))
-	    {	
+	    {
 		x=copysign(x,one);
 		y=copysign(y,one);
-		if(y > x) 
+		if(y > x)
 		    { temp=x; x=y; y=temp; }
 		if(x == zero) return(zero);
 		if(y == zero) return(x);
 		exp= logb(x);
 		x=scalb(x,-exp);
-		if(exp-(int)logb(y) > ibig ) 
+		if(exp-(int)logb(y) > ibig )
 			/* raise inexact flag and return |x| */
 		   { one+small; return(scalb(x,exp)); }
 		else y=scalb(y,-exp);
@@ -217,7 +217,7 @@ double x, y;
 
 	    else if(y==y)   	   /* y is +-INF */
 		     return(copysign(y,one));
-	    else 
+	    else
 		     return(y);	   /* y is NaN and x is finite */
 
 	else if(x==x) 		   /* x is +-INF */
diff --git a/lib/libm/ieee/cbrt.c b/lib/libm/ieee/cbrt.c
index fe5fb95..6d73779 100644
--- a/lib/libm/ieee/cbrt.c
+++ b/lib/libm/ieee/cbrt.c
@@ -50,8 +50,8 @@ static char sccsid[] = "@(#)cbrt.c	8.1 (Berkeley) 6/4/93";
  * long interger at the address of a floating point number will be the
  * leading 32 bits of that floating point number (i.e., sign, exponent,
  * and the 20 most significant bits).
- * On a National machine, it has different ordering; therefore, this code 
- * must be compiled with flag -DNATIONAL. 
+ * On a National machine, it has different ordering; therefore, this code
+ * must be compiled with flag -DNATIONAL.
  */
 #if !defined(vax)&&!defined(tahoe)
 
@@ -65,7 +65,7 @@ static const double
 	    F= 45./28.,
 	    G= 5./14.;
 
-double cbrt(x) 
+double cbrt(x)
 double x;
 {
 	double r,s,t=0.0,w;
@@ -93,15 +93,15 @@ double x;
 	   t*=x; pt[n0]=pt[n0]/3+B2;
 	  }
 	else
-	  pt[n0]=px[n0]/3+B1;	
+	  pt[n0]=px[n0]/3+B1;
 
 
     /* new cbrt to 23 bits, may be implemented in single precision */
 	r=t*t/x;
 	s=C+r*t;
-	t*=G+F/(s+E+D/s);	
+	t*=G+F/(s+E+D/s);
 
-    /* chopped to 20 bits and make it larger than cbrt(x) */ 
+    /* chopped to 20 bits and make it larger than cbrt(x) */
 	pt[n1]=0; pt[n0]+=0x00000001;
 
 
diff --git a/lib/libm/ieee/support.c b/lib/libm/ieee/support.c
index e976839..0ae10bb 100644
--- a/lib/libm/ieee/support.c
+++ b/lib/libm/ieee/support.c
@@ -35,37 +35,37 @@
 static char sccsid[] = "@(#)support.c	8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
-/* 
- * Some IEEE standard 754 recommended functions and remainder and sqrt for 
+/*
+ * Some IEEE standard 754 recommended functions and remainder and sqrt for
  * supporting the C elementary functions.
  ******************************************************************************
  * WARNING:
  *      These codes are developed (in double) to support the C elementary
  * functions temporarily. They are not universal, and some of them are very
- * slow (in particular, drem and sqrt is extremely inefficient). Each 
- * computer system should have its implementation of these functions using 
+ * slow (in particular, drem and sqrt is extremely inefficient). Each
+ * computer system should have its implementation of these functions using
  * its own assembler.
  ******************************************************************************
  *
  * IEEE 754 required operations:
- *     drem(x,p) 
+ *     drem(x,p)
  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
  *              nearest x/y; in half way case, choose the even one.
- *     sqrt(x) 
- *              returns the square root of x correctly rounded according to 
+ *     sqrt(x)
+ *              returns the square root of x correctly rounded according to
  *		the rounding mod.
  *
  * IEEE 754 recommended functions:
- * (a) copysign(x,y) 
- *              returns x with the sign of y. 
- * (b) scalb(x,N) 
+ * (a) copysign(x,y)
+ *              returns x with the sign of y.
+ * (b) scalb(x,N)
  *              returns  x * (2**N), for integer values N.
- * (c) logb(x) 
- *              returns the unbiased exponent of x, a signed integer in 
- *              double precision, except that logb(0) is -INF, logb(INF) 
+ * (c) logb(x)
+ *              returns the unbiased exponent of x, a signed integer in
+ *              double precision, except that logb(0) is -INF, logb(INF)
  *              is +INF, and logb(NAN) is that NAN.
- * (d) finite(x) 
- *              returns the value TRUE if -INF < x < +INF and returns 
+ * (d) finite(x)
+ *              returns the value TRUE if -INF < x < +INF and returns
  *              FALSE otherwise.
  *
  *
@@ -78,7 +78,7 @@ static char sccsid[] = "@(#)support.c	8.1 (Berkeley) 6/4/93";
 #if defined(vax)||defined(tahoe)      /* VAX D format */
 #include <errno.h>
     static const unsigned short msign=0x7fff , mexp =0x7f80 ;
-    static const short  prep1=57, gap=7, bias=129           ;   
+    static const short  prep1=57, gap=7, bias=129           ;
     static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
 #else	/* defined(vax)||defined(tahoe) */
     static const unsigned short msign=0x7fff, mexp =0x7ff0  ;
@@ -97,7 +97,7 @@ double x; int N;
         unsigned short *px=(unsigned short *) &x;
 #endif	/* national */
 
-        if( x == zero )  return(x); 
+        if( x == zero )  return(x);
 
 #if defined(vax)||defined(tahoe)
         if( (k= *px & mexp ) != ~msign ) {
@@ -117,7 +117,7 @@ double x; int N;
                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
                 else x=novf+novf;               /* overflow */
             else
-                if( k > -prep1 ) 
+                if( k > -prep1 )
                                         /* gradual underflow */
                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
                 else
@@ -147,7 +147,7 @@ double x,y;
 }
 
 double logb(x)
-double x; 
+double x;
 {
 
 #ifdef national
@@ -164,8 +164,8 @@ double x;
                 return ( (k>>gap) - bias );
             else if( x != zero)
                 return ( -1022.0 );
-            else        
-                return(-(1.0/zero));    
+            else
+                return(-(1.0/zero));
         else if(x != x)
             return(x);
         else
@@ -174,7 +174,7 @@ double x;
 }
 
 finite(x)
-double x;    
+double x;
 {
 #if defined(vax)||defined(tahoe)
         return(1);
@@ -192,16 +192,16 @@ double x,p;
 {
         short sign;
         double hp,dp,tmp;
-        unsigned short  k; 
+        unsigned short  k;
 #ifdef national
         unsigned short
-              *px=(unsigned short *) &x  +3, 
+              *px=(unsigned short *) &x  +3,
               *pp=(unsigned short *) &p  +3,
               *pd=(unsigned short *) &dp +3,
               *pt=(unsigned short *) &tmp+3;
 #else	/* national */
         unsigned short
-              *px=(unsigned short *) &x  , 
+              *px=(unsigned short *) &x  ,
               *pp=(unsigned short *) &p  ,
               *pd=(unsigned short *) &dp ,
               *pt=(unsigned short *) &tmp;
@@ -230,13 +230,13 @@ double x,p;
 #endif	/* defined(vax)||defined(tahoe) */
 		{ if (p != p) return p; else return x;}
 
-        else  if ( ((*pp & mexp)>>gap) <= 1 ) 
+        else  if ( ((*pp & mexp)>>gap) <= 1 )
                 /* subnormal p, or almost subnormal p */
             { double b; b=scalb(1.0,(int)prep1);
               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
         else  if ( p >= novf/2)
             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
-        else 
+        else
             {
                 dp=p+p; hp=p/2;
                 sign= *px & ~msign ;
@@ -294,13 +294,13 @@ double x;
 	}
 
     /* sqrt(INF) is INF */
-        if(!finite(x)) return(x);               
+        if(!finite(x)) return(x);
 
     /* scale x to [1,4) */
         n=logb(x);
         x=scalb(x,-n);
         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
-        m += n; 
+        m += n;
         n = m/2;
         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
 
@@ -313,23 +313,23 @@ double x;
                 else
                     s *= 2;
                 }
-            
+
     /* generate the last bit and determine the final rounding */
-            r/=2; x *= 4; 
+            r/=2; x *= 4;
             if(x==zero) goto end; 100+r; /* trigger inexact flag */
             if(s<x) {
                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
-                t = (x-s)-5; 
+                t = (x-s)-5;
                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
                 b=1.0+r/4;   if(b>1.0) t=1;	/* b>1 : Round-to-(+INF) */
                 if(t>=0) q+=r; }	      /* else: Round-to-nearest */
-            else { 
-                s *= 2; x *= 4; 
-                t = (x-s)-1; 
+            else {
+                s *= 2; x *= 4;
+                t = (x-s)-1;
                 b=1.0+3*r/4; if(b==1.0) goto end;
                 b=1.0+r/4;   if(b>1.0) t=1;
                 if(t>=0) q+=r; }
-            
+
 end:        return(scalb(q,n));
 }
 
@@ -342,13 +342,13 @@ end:        return(scalb(q,n));
  *
  * Warning: this code should not get compiled in unless ALL of
  * the following machine-dependent routines are supplied.
- * 
+ *
  * Required machine dependent functions (not on a VAX):
  *     swapINX(i): save inexact flag and reset it to "i"
  *     swapENI(e): save inexact enable and reset it to "e"
  */
 
-double drem(x,y)	
+double drem(x,y)
 double x,y;
 {
 
@@ -363,8 +363,8 @@ double x,y;
 	double hy,y1,t,t1;
 	short k;
 	long n;
-	int i,e; 
-	unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
+	int i,e;
+	unsigned short xexp,yexp, *px  =(unsigned short *) &x  ,
 	      		nx,nf,	  *py  =(unsigned short *) &y  ,
 	      		sign,	  *pt  =(unsigned short *) &t  ,
 	      			  *pt1 =(unsigned short *) &t1 ;
@@ -381,7 +381,7 @@ double x,y;
 /* save the inexact flag and inexact enable in i and e respectively
  * and reset them to zero
  */
-	i=swapINX(0);	e=swapENI(0);	
+	i=swapINX(0);	e=swapENI(0);
 
 /* subnormal number */
 	nx=0;
@@ -391,7 +391,7 @@ double x,y;
 	if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
 
 	nf=nx;
-	py[n0] &= 0x7fff;	
+	py[n0] &= 0x7fff;
 	px[n0] &= 0x7fff;
 
 /* mask off the least significant 27 bits of y */
@@ -408,7 +408,7 @@ loop:
 	    if(k>0) 	/* if x/y >= 2**26, scale up y so that x/y < 2**26 */
 		{pt[n0]+=k;pt1[n0]+=k;}
 	    n=x/t; x=(x-n*t1)-n*(t-t1);
-	}	
+	}
     /* end while (x > y) */
 
 	if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
@@ -416,14 +416,14 @@ loop:
 /* final adjustment */
 
 	hy=y/2.0;
-	if(x>hy||((x==hy)&&n%2==1)) x-=y; 
+	if(x>hy||((x==hy)&&n%2==1)) x-=y;
 	px[n0] ^= sign;
 	if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
 
 /* restore inexact flag and inexact enable */
-	swapINX(i); swapENI(e);	
+	swapINX(i); swapENI(e);
 
-	return(x);	
+	return(x);
 }
 #endif
 
@@ -435,7 +435,7 @@ loop:
  *
  * Warning: this code should not get compiled in unless ALL of
  * the following machine-dependent routines are supplied.
- * 
+ *
  * Required machine dependent functions:
  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
@@ -456,16 +456,16 @@ double x;
 	double const b54=134217728.*134217728.; /* b54=2**54 */
         long mx,scalx;
 	long const mexp=0x7ff00000;
-        int i,j,r,e,swapINX(),swapRM(),swapENI();       
+        int i,j,r,e,swapINX(),swapRM(),swapENI();
         unsigned long *py=(unsigned long *) &y   ,
                       *pt=(unsigned long *) &t   ,
                       *px=(unsigned long *) &x   ;
 #ifdef national         /* ordering of word in a floating point number */
-        const int n0=1, n1=0; 
+        const int n0=1, n1=0;
 #else
-        const int n0=0, n1=1; 
+        const int n0=0, n1=1;
 #endif
-/* Rounding Mode:  RN ...round-to-nearest 
+/* Rounding Mode:  RN ...round-to-nearest
  *                 RZ ...round-towards 0
  *                 RP ...round-towards +INF
  *		   RM ...round-towards -INF
@@ -502,10 +502,10 @@ double x;
         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
 
 /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
-        t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
+        t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
 
-/* twiddle last bit to force y correctly rounded */ 
+/* twiddle last bit to force y correctly rounded */
         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
         swapINX(0);     /* ...clear INEXACT flag */
         swapENI(e);     /* ...restore inexact enable status */