Trim off whitespace

libm/Makefile

@@ -6,14 +6,14 @@
 # math library for Apple's MacOS X/Darwin math library, which was
 # itself swiped from FreeBSD.  The original copyright information
 # is as follows:
-# 
+#
 #     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.
-# 
+#
 # It has been ported to work with uClibc and generally behave
 # by Erik Andersen <andersen@codepoet.org>
 #
@@ -36,7 +36,7 @@ include $(TOPDIR)Rules.mak
 
 CFLAGS+=$(SSP_ALL_CFLAGS)
 
-DIRS = 
+DIRS =
 ifeq ($(strip $(HAS_FPU)),y)
 ifeq ($(TARGET_ARCH),$(wildcard $(TARGET_ARCH)))
 DIRS = $(TARGET_ARCH)
@@ -71,7 +71,7 @@ else
 CSRC =   w_acos.c w_asin.c s_atan.c w_atan2.c s_ceil.c s_cos.c \
 	 w_cosh.c w_exp.c s_fabs.c s_floor.c w_fmod.c s_frexp.c \
 	 s_ldexp.c w_log.c w_log10.c s_modf.c w_pow.c s_sin.c \
-	 w_sinh.c w_sqrt.c s_tan.c s_tanh.c 
+	 w_sinh.c w_sqrt.c s_tan.c s_tanh.c
 CSRC+=   s_expm1.c s_scalbn.c s_copysign.c e_acos.c e_asin.c e_atan2.c \
 	 k_cos.c e_cosh.c e_exp.c e_fmod.c e_log.c e_log10.c e_pow.c \
 	 k_sin.c e_sinh.c e_sqrt.c k_tan.c e_rem_pio2.c k_rem_pio2.c \

libm/README

@@ -7,10 +7,10 @@ is as follows:
 
 	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.
 
-It has been ported to work with uClibc and generally behave 
+It has been ported to work with uClibc and generally behave
 by Erik Andersen <andersen@codepoet.org>
   22 May, 2001
 

libm/e_acos.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,14 +15,14 @@ static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $";
 #endif
 
 /* __ieee754_acos(x)
- * Method :                  
+ * Method :
  *	acos(x)  = pi/2 - asin(x)
  *	acos(-x) = pi/2 + asin(x)
  * For |x|<=0.5
  *	acos(x) = pi/2 - (x + x*x^2*R(x^2))	(see asin.c)
  * For x>0.5
  * 	acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
- *		= 2asin(sqrt((1-x)/2))  
+ *		= 2asin(sqrt((1-x)/2))
  *		= 2s + 2s*z*R(z) 	...z=(1-x)/2, s=sqrt(z)
  *		= 2f + (2c + 2s*z*R(z))
  *     where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
@@ -42,9 +42,9 @@ static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $";
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 one=  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 pi =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */

libm/e_acosh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -16,7 +16,7 @@ static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $";
 
 /* __ieee754_acosh(x)
  * Method :
- *	Based on 
+ *	Based on
  *		acosh(x) = log [ x + sqrt(x*x-1) ]
  *	we have
  *		acosh(x) := log(x)+ln2,	if x is large; else
@@ -32,9 +32,9 @@ static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $";
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 one	= 1.0,
 ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
@@ -45,7 +45,7 @@ ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
 	double __ieee754_acosh(x)
 	double x;
 #endif
-{	
+{
 	double t;
 	int32_t hx;
 	u_int32_t lx;
@@ -55,7 +55,7 @@ ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
 	} else if(hx >=0x41b00000) {	/* x > 2**28 */
 	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
 	        return x+x;
-	    } else 
+	    } else
 		return __ieee754_log(x)+ln2;	/* acosh(huge)=log(2x) */
 	} else if(((hx-0x3ff00000)|lx)==0) {
 	    return 0.0;			/* acosh(1) = 0 */

libm/e_asin.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,12 +15,12 @@ static char rcsid[] = "$NetBSD: e_asin.c,v 1.9 1995/05/12 04:57:22 jtc Exp $";
 #endif
 
 /* __ieee754_asin(x)
- * Method :                  
+ * Method :
  *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
  *	we approximate asin(x) on [0,0.5] by
  *		asin(x) = x + x*x^2*R(x^2)
  *	where
- *		R(x^2) is a rational approximation of (asin(x)-x)/x^3 
+ *		R(x^2) is a rational approximation of (asin(x)-x)/x^3
  *	and its remez error is bounded by
  *		|(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
  *
@@ -49,9 +49,9 @@ static char rcsid[] = "$NetBSD: e_asin.c,v 1.9 1995/05/12 04:57:22 jtc Exp $";
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 huge =  1.000e+300,
@@ -86,12 +86,12 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 	    GET_LOW_WORD(lx,x);
 	    if(((ix-0x3ff00000)|lx)==0)
 		    /* asin(1)=+-pi/2 with inexact */
-		return x*pio2_hi+x*pio2_lo;	
-	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
+		return x*pio2_hi+x*pio2_lo;
+	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */
 	} else if (ix<0x3fe00000) {	/* |x|<0.5 */
 	    if(ix<0x3e400000) {		/* if |x| < 2**-27 */
 		if(huge+x>one) return x;/* return x with inexact if x!=0*/
-	    } else { 
+	    } else {
 		t = x*x;
 		p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
 		q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
@@ -116,6 +116,6 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 	    p  = 2.0*s*r-(pio2_lo-2.0*c);
 	    q  = pio4_hi-2.0*w;
 	    t  = pio4_hi-(p-q);
-	}    
-	if(hx>0) return t; else return -t;    
+	}
+	if(hx>0) return t; else return -t;
 }

libm/e_atan2.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -17,7 +17,7 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $";
 /* __ieee754_atan2(y,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,
  *
@@ -35,9 +35,9 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $";
  *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
- * compiler will convert from decimal to binary accurately enough 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
 
@@ -45,9 +45,9 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $";
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 tiny  = 1.0e-300,
 zero  = 0.0,
@@ -62,7 +62,7 @@ pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
 	double __ieee754_atan2(y,x)
 	double  y,x;
 #endif
-{  
+{
 	double z;
 	int32_t k,m,hx,hy,ix,iy;
 	u_int32_t lx,ly;
@@ -80,7 +80,7 @@ pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
     /* when y = 0 */
 	if((iy|ly)==0) {
 	    switch(m) {
-		case 0: 
+		case 0:
 		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
 		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
 		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
@@ -88,7 +88,7 @@ pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
 	}
     /* when x = 0 */
 	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
-	    
+
     /* when x is INF */
 	if(ix==0x7ff00000) {
 	    if(iy==0x7ff00000) {

libm/e_atanh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -21,7 +21,7 @@ static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $";
  *                  1              2x                          x
  *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
  *                  2             1 - x                      1 - x
- *	
+ *
  * 	For x<0.5
  *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
  *
@@ -61,14 +61,14 @@ static double zero = 0.0;
 	ix = hx&0x7fffffff;
 	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
 	    return (x-x)/(x-x);
-	if(ix==0x3ff00000) 
+	if(ix==0x3ff00000)
 	    return x/zero;
 	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
 	SET_HIGH_WORD(x,ix);
 	if(ix<0x3fe00000) {		/* x < 0.5 */
 	    t = x+x;
 	    t = 0.5*log1p(t+t*x/(one-x));
-	} else 
+	} else
 	    t = 0.5*log1p((x+x)/(one-x));
 	if(hx>=0) return t; else return -t;
 }

libm/e_cosh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,18 +15,18 @@ static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $";
 #endif
 
 /* __ieee754_cosh(x)
- * Method : 
+ * Method :
  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
- *	1. Replace x by |x| (cosh(x) = cosh(-x)). 
- *	2. 
- *		                                        [ exp(x) - 1 ]^2 
+ *	1. Replace x by |x| (cosh(x) = cosh(-x)).
+ *	2.
+ *		                                        [ exp(x) - 1 ]^2
  *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
  *			       			           2*exp(x)
  *
  *		                                  exp(x) +  1/exp(x)
  *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
  *			       			          2
- *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2 
+ *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2
  *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
  *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
  *
@@ -50,7 +50,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
 	double __ieee754_cosh(x)
 	double x;
 #endif
-{	
+{
 	double t,w;
 	int32_t ix;
 	u_int32_t lx;
@@ -60,7 +60,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
 	ix &= 0x7fffffff;
 
     /* x is INF or NaN */
-	if(ix>=0x7ff00000) return x*x;	
+	if(ix>=0x7ff00000) return x*x;
 
     /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
 	if(ix<0x3fd62e43) {
@@ -81,7 +81,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300;
 
     /* |x| in [log(maxdouble), overflowthresold] */
 	GET_LOW_WORD(lx,x);
-	if (ix<0x408633CE || 
+	if (ix<0x408633CE ||
 	      ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
 	    w = __ieee754_exp(half*fabs(x));
 	    t = half*w;

libm/e_exp.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -22,36 +22,36 @@ static char rcsid[] = "$NetBSD: e_exp.c,v 1.8 1995/05/10 20:45:03 jtc Exp $";
  *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
  *	Given x, find r and integer k such that
  *
- *               x = k*ln2 + r,  |r| <= 0.5*ln2.  
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2.
  *
- *      Here r will be represented as r = hi-lo for better 
+ *      Here r will be represented as r = hi-lo for better
  *	accuracy.
  *
  *   2. Approximation of exp(r) by a special rational function on
  *	the interval [0,0.34658]:
  *	Write
  *	    R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
- *      We use a special Reme algorithm on [0,0.34658] to generate 
- * 	a polynomial of degree 5 to approximate R. The maximum error 
+ *      We use a special Reme algorithm on [0,0.34658] to generate
+ * 	a polynomial of degree 5 to approximate R. The maximum error
  *	of this polynomial approximation is bounded by 2**-59. In
  *	other words,
  *	    R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
  *  	(where z=r*r, and the values of P1 to P5 are listed below)
  *	and
  *	    |                  5          |     -59
- *	    | 2.0+P1*z+...+P5*z   -  R(z) | <= 2 
+ *	    | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
  *	    |                             |
  *	The computation of exp(r) thus becomes
  *                             2*r
  *		exp(r) = 1 + -------
  *		              R - r
- *                                 r*R1(r)	
+ *                                 r*R1(r)
  *		       = 1 + r + ----------- (for better accuracy)
  *		                  2 - R1(r)
  *	where
  *			         2       4             10
  *		R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
- *	
+ *
  *   3. Scale back to obtain exp(x):
  *	From step 1, we have
  *	   exp(x) = 2^k * exp(r)
@@ -66,13 +66,13 @@ static char rcsid[] = "$NetBSD: e_exp.c,v 1.8 1995/05/10 20:45:03 jtc Exp $";
  *	1 ulp (unit in the last place).
  *
  * Misc. info.
- *	For IEEE double 
+ *	For IEEE double
  *	    if x >  7.09782712893383973096e+02 then exp(x) overflow
  *	    if x < -7.45133219101941108420e+02 then exp(x) underflow
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
  * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
@@ -128,7 +128,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
             if(hx>=0x7ff00000) {
 	        u_int32_t lx;
 		GET_LOW_WORD(lx,x);
-		if(((hx&0xfffff)|lx)!=0) 
+		if(((hx&0xfffff)|lx)!=0)
 		     return x+x; 		/* NaN */
 		else return (xsb==0)? x:0.0;	/* exp(+-inf)={inf,0} */
 	    }
@@ -137,7 +137,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 	}
 
     /* argument reduction */
-	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */ 
+	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */
 	    if(hx < 0x3FF0A2B2) {	/* and |x| < 1.5 ln2 */
 		hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
 	    } else {
@@ -147,7 +147,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
 		lo = t*ln2LO[0];
 	    }
 	    x  = hi - lo;
-	} 
+	}
 	else if(hx < 0x3e300000)  {	/* when |x|<2**-28 */
 	    if(huge+x>one) return one+x;/* trigger inexact */
 	}
@@ -156,7 +156,7 @@ P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
     /* x is now in primary range */
 	t  = x*x;
 	c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
-	if(k==0) 	return one-((x*c)/(c-2.0)-x); 
+	if(k==0) 	return one-((x*c)/(c-2.0)-x);
 	else 		y = one-((lo-(x*c)/(2.0-c))-hi);
 	if(k >= -1021) {
 	    u_int32_t hy;

libm/e_fmod.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $";
 #endif
 
-/* 
+/*
  * __ieee754_fmod(x,y)
  * Return x mod y in exact arithmetic
  * Method: shift and subtract
@@ -51,7 +51,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 	    return (x*y)/(x*y);
 	if(hx<=hy) {
 	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
-	    if(lx==ly) 
+	    if(lx==ly)
 		return Zero[(u_int32_t)sx>>31];	/* |x|=|y| return x*0*/
 	}
 
@@ -74,7 +74,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 	} else iy = (hy>>20)-1023;
 
     /* set up {hx,lx}, {hy,ly} and align y to x */
-	if(ix >= -1022) 
+	if(ix >= -1022)
 	    hx = 0x00100000|(0x000fffff&hx);
 	else {		/* subnormal x, shift x to normal */
 	    n = -1022-ix;
@@ -86,7 +86,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 		lx = 0;
 	    }
 	}
-	if(iy >= -1022) 
+	if(iy >= -1022)
 	    hy = 0x00100000|(0x000fffff&hy);
 	else {		/* subnormal y, shift y to normal */
 	    n = -1022-iy;
@@ -115,7 +115,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,};
 
     /* convert back to floating value and restore the sign */
 	if((hx|lx)==0) 			/* return sign(x)*0 */
-	    return Zero[(u_int32_t)sx>>31];	
+	    return Zero[(u_int32_t)sx>>31];
 	while(hx<0x00100000) {		/* normalize x */
 	    hx = hx+hx+(lx>>31); lx = lx+lx;
 	    iy -= 1;

libm/e_gamma.c

@@ -6,7 +6,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  *

libm/e_gamma_r.c

@@ -6,15 +6,15 @@
  *
  * 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.
  * ====================================================
  *
  */
 
 /* __ieee754_gamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
  *
  * Method: See __ieee754_lgamma_r
  */

libm/e_hypot.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -16,12 +16,12 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
 
 /* __ieee754_hypot(x,y)
  *
- * Method :                  
- *	If (assume round-to-nearest) z=x*x+y*y 
- *	has error less than sqrt(2)/2 ulp, than 
+ * Method :
+ *	If (assume round-to-nearest) z=x*x+y*y
+ *	has error less than sqrt(2)/2 ulp, than
  *	sqrt(z) has error less than 1 ulp (exercise).
  *
- *	So, compute sqrt(x*x+y*y) with some care as 
+ *	So, compute sqrt(x*x+y*y) with some care as
  *	follows to get the error below 1 ulp:
  *
  *	Assume x>y>0;
@@ -31,10 +31,10 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
  *	where x1 = x with lower 32 bits cleared, x2 = x-x1; else
  *	2. if x <= 2y use
  *		t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
- *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 
+ *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
  *	y1= y with lower 32 bits chopped, y2 = y-y1.
- *		
- *	NOTE: scaling may be necessary if some argument is too 
+ *
+ *	NOTE: scaling may be necessary if some argument is too
  *	      large or too tiny
  *
  * Special cases:
@@ -42,8 +42,8 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
  *	hypot(x,y) is NAN if x or y is NAN.
  *
  * Accuracy:
- * 	hypot(x,y) returns sqrt(x^2+y^2) with error less 
- * 	than 1 ulps (units in the last place) 
+ * 	hypot(x,y) returns sqrt(x^2+y^2) with error less
+ * 	than 1 ulps (units in the last place)
  */
 
 #include "math.h"
@@ -84,7 +84,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $";
 	   SET_HIGH_WORD(b,hb);
 	}
 	if(hb < 0x20b00000) {	/* b < 2**-500 */
-	    if(hb <= 0x000fffff) {	/* subnormal b or 0 */	
+	    if(hb <= 0x000fffff) {	/* subnormal b or 0 */
 	        u_int32_t low;
 		GET_LOW_WORD(low,b);
 		if((hb|low)==0) return a;

libm/e_j0.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -33,20 +33,20 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $";
  * 	   (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) = u00 + u01*z + ... + u06*z^6
  *		V(z) = 1  + v01*z + ... + v04*z^4
  *	   with absolute approximation error bounded by 2**-72.
@@ -69,9 +69,9 @@ static double pzero(), qzero();
 #endif
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 huge 	= 1e300,
 one	= 1.0,
@@ -94,9 +94,9 @@ static double zero = 0.0;
 #endif
 
 #ifdef __STDC__
-	double __ieee754_j0(double x) 
+	double __ieee754_j0(double x)
 #else
-	double __ieee754_j0(x) 
+	double __ieee754_j0(x)
 	double x;
 #endif
 {
@@ -163,9 +163,9 @@ v03  =  2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */
 v04  =  4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
 
 #ifdef __STDC__
-	double __ieee754_y0(double x) 
+	double __ieee754_y0(double x)
 #else
-	double __ieee754_y0(x) 
+	double __ieee754_y0(x)
 	double x;
 #endif
 {
@@ -175,7 +175,7 @@ v04  =  4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
     /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7ff00000) return  one/(x+x*x); 
+	if(ix>=0x7ff00000) return  one/(x+x*x);
         if((ix|lx)==0) return -one/zero;
         if(hx<0) return zero/zero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
@@ -349,7 +349,7 @@ static double pS2[5] = {
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
 	return one+ r/s;
 }
-		
+
 
 /* For x >= 8, the asymptotic expansions of qzero is
  *	-1/8 s + 75/1024 s^3 - ..., where s = 1/x.

libm/e_j1.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -34,16 +34,16 @@ static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $";
  * 	   (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,
@@ -69,9 +69,9 @@ static double pone(), qone();
 #endif
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 huge    = 1e300,
 one	= 1.0,
@@ -95,9 +95,9 @@ static double zero    = 0.0;
 #endif
 
 #ifdef __STDC__
-	double __ieee754_j1(double x) 
+	double __ieee754_j1(double x)
 #else
-	double __ieee754_j1(x) 
+	double __ieee754_j1(x)
 	double x;
 #endif
 {
@@ -164,9 +164,9 @@ static double V0[5] = {
 };
 
 #ifdef __STDC__
-	double __ieee754_y1(double x) 
+	double __ieee754_y1(double x)
 #else
-	double __ieee754_y1(x) 
+	double __ieee754_y1(x)
 	double x;
 #endif
 {
@@ -176,7 +176,7 @@ static double V0[5] = {
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
     /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7ff00000) return  one/(x+x*x); 
+	if(ix>=0x7ff00000) return  one/(x+x*x);
         if((ix|lx)==0) return -one/zero;
         if(hx<0) return zero/zero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
@@ -206,10 +206,10 @@ static double V0[5] = {
                     z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
                 }
                 return z;
-        } 
+        }
         if(ix<=0x3c900000) {    /* 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]))));
@@ -347,7 +347,7 @@ static double ps2[5] = {
         s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
         return one+ r/s;
 }
-		
+
 
 /* For x >= 8, the asymptotic expansions of qone is
  *	3/8 s - 105/1024 s^3 - ..., where s = 1/x.

libm/e_jn.c

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

libm/e_lgamma.c

@@ -6,7 +6,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  *

libm/e_lgamma_r.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,12 +15,12 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $
 #endif
 
 /* __ieee754_lgamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
  *
  * Method:
  *   1. Argument Reduction for 0 < x <= 8
- * 	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may 
+ * 	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
  * 	reduce x to a number in [1.5,2.5] by
  * 		lgamma(1+s) = log(s) + lgamma(s)
  *	for example,
@@ -58,36 +58,36 @@ static char rcsid[] = "$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $
  *	by
  *	  			    3       5             11
  *		w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
- *	where 
+ *	where
  *		|w - f(z)| < 2**-58.74
- *		
+ *
  *   4. For negative x, since (G is gamma function)
  *		-x*G(-x)*G(x) = pi/sin(pi*x),
  * 	we have
  * 		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
  *	since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
- *	Hence, for x<0, signgam = sign(sin(pi*x)) and 
+ *	Hence, for x<0, signgam = sign(sin(pi*x)) and
  *		lgamma(x) = log(|Gamma(x)|)
  *			  = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
- *	Note: one should avoid compute pi*(-x) directly in the 
+ *	Note: one should avoid compute pi*(-x) directly in the
  *	      computation of sin(pi*(-x)).
- *		
+ *
  *   5. Special Cases
  *		lgamma(2+s) ~ s*(1-Euler) for tiny s
  *		lgamma(1)=lgamma(2)=0
  *		lgamma(x) ~ -log(x) for tiny x
  *		lgamma(0) = lgamma(inf) = inf
  *	 	lgamma(-integer) = +-inf
- *	
+ *
  */
 
 #include "math.h"
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 two52=  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
@@ -204,9 +204,9 @@ __inline__
         }
 	switch (n) {
 	    case 0:   y =  __kernel_sin(pi*y,zero,0); break;
-	    case 1:   
+	    case 1:
 	    case 2:   y =  __kernel_cos(pi*(0.5-y),zero); break;
-	    case 3:  
+	    case 3:
 	    case 4:   y =  __kernel_sin(pi*(one-y),zero,0); break;
 	    case 5:
 	    case 6:   y = -__kernel_cos(pi*(y-1.5),zero); break;
@@ -279,7 +279,7 @@ __inline__
 		p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
 		p  = z*p1-(tt-w*(p2+y*p3));
 		r += (tf + p); break;
-	      case 2:	
+	      case 2:
 		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
 		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
 		r += (-0.5*y + p1/p2);
@@ -308,7 +308,7 @@ __inline__
 	    y = z*z;
 	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
 	    r = (x-half)*(t-one)+w;
-	} else 
+	} else
     /* 2**58 <= x <= inf */
 	    r =  x*(__ieee754_log(x)-one);
 	if(hx<0) r = nadj - r;

libm/e_log.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -17,17 +17,17 @@ static char rcsid[] = "$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $";
 /* __ieee754_log(x)
  * Return the logrithm of x
  *
- * Method :                  
- *   1. Argument Reduction: find k and f such that 
- *			x = 2^k * (1+f), 
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *			x = 2^k * (1+f),
  *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
  *
  *   2. Approximation of log(1+f).
  *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  *	     	 = 2s + s*R
- *      We use a special Reme algorithm on [0,0.1716] to generate 
- * 	a polynomial of degree 14 to approximate R The maximum error 
+ *      We use a special Reme algorithm on [0,0.1716] to generate
+ * 	a polynomial of degree 14 to approximate R The maximum error
  *	of this polynomial approximation is bounded by 2**-58.45. In
  *	other words,
  *		        2      4      6      8      10      12      14
@@ -35,22 +35,22 @@ static char rcsid[] = "$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $";
  *  	(the values of Lg1 to Lg7 are listed in the program)
  *	and
  *	    |      2          14          |     -58.45
- *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
+ *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2
  *	    |                             |
  *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
  *	In order to guarantee error in log below 1ulp, we compute log
  *	by
  *		log(1+f) = f - s*(f - R)	(if f is not too large)
  *		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
- *	
- *	3. Finally,  log(x) = k*ln2 + log(1+f).  
+ *
+ *	3. Finally,  log(x) = k*ln2 + log(1+f).
  *			    = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
- *	   Here ln2 is split into two floating point number: 
+ *	   Here ln2 is split into two floating point number:
  *			ln2_hi + ln2_lo,
  *	   where n*ln2_hi is always exact for |n| < 2000.
  *
  * Special cases:
- *	log(x) is NaN with signal if x < 0 (including -INF) ; 
+ *	log(x) is NaN with signal if x < 0 (including -INF) ;
  *	log(+INF) is +INF; log(0) is -INF with signal;
  *	log(NaN) is that NaN with no signal.
  *
@@ -59,9 +59,9 @@ static char rcsid[] = "$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $";
  *	1 ulp (unit in the last place).
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
- * compiler will convert from decimal to binary accurately enough 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
 
@@ -105,12 +105,12 @@ static double zero   =  0.0;
 
 	k=0;
 	if (hx < 0x00100000) {			/* x < 2**-1022  */
-	    if (((hx&0x7fffffff)|lx)==0) 
+	    if (((hx&0x7fffffff)|lx)==0)
 		return -two54/zero;		/* log(+-0)=-inf */
 	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
 	    k -= 54; x *= two54; /* subnormal number, scale up x */
 	    GET_HIGH_WORD(hx,x);
-	} 
+	}
 	if (hx >= 0x7ff00000) return x+x;
 	k += (hx>>20)-1023;
 	hx &= 0x000fffff;
@@ -126,14 +126,14 @@ static double zero   =  0.0;
 	    if(k==0) return f-R; else {dk=(double)k;
 	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
 	}
- 	s = f/(2.0+f); 
+ 	s = f/(2.0+f);
 	dk = (double)k;
 	z = s*s;
 	i = hx-0x6147a;
 	w = z*z;
 	j = 0x6b851-hx;
-	t1= w*(Lg2+w*(Lg4+w*Lg6)); 
-	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 
+	t1= w*(Lg2+w*(Lg4+w*Lg6));
+	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
 	i |= j;
 	R = t2+t1;
 	if(i>0) {

libm/e_log10.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -16,26 +16,26 @@ static char rcsid[] = "$NetBSD: e_log10.c,v 1.9 1995/05/10 20:45:51 jtc Exp $";
 
 /* __ieee754_log10(x)
  * Return the base 10 logarithm of x
- * 
+ *
  * Method :
  *	Let log10_2hi = leading 40 bits of log10(2) and
  *	    log10_2lo = log10(2) - log10_2hi,
  *	    ivln10   = 1/log(10) rounded.
  *	Then
- *		n = ilogb(x), 
+ *		n = ilogb(x),
  *		if(n<0)  n = n+1;
  *		x = scalbn(x,-n);
  *		log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
  *
  * Note 1:
- *	To guarantee log10(10**n)=n, where 10**n is normal, the rounding 
+ *	To guarantee log10(10**n)=n, where 10**n is normal, the rounding
  *	mode must set to Round-to-Nearest.
  * Note 2:
  *	[1/log(10)] rounded to 53 bits has error  .198   ulps;
  *	log10 is monotonic at all binary break points.
  *
  * 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;
  *	log10(10**N) = N  for N=0,1,...,22.

libm/e_pow.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -21,7 +21,7 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
  *	1. Compute and return log2(x) in two pieces:
  *		log2(x) = w1 + w2,
  *	   where w1 has 53-24 = 29 bit trailing zeros.
- *	2. Perform y*log2(x) = n+y' by simulating muti-precision 
+ *	2. Perform y*log2(x) = n+y' by simulating muti-precision
  *	   arithmetic, where |y'|<=0.5.
  *	3. Return x**y = 2**n*exp(y'*log2)
  *
@@ -49,13 +49,13 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
  * Accuracy:
  *	pow(x,y) returns x**y nearly rounded. In particular
  *			pow(integer,integer)
- *	always returns the correct integer provided it is 
+ *	always returns the correct integer provided it is
  *	representable.
  *
  * Constants :
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
- * compiler will convert from decimal to binary accurately enough 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
 
@@ -63,9 +63,9 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 bp[] = {1.0, 1.5,},
 dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
@@ -117,12 +117,12 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
 
     /* y==zero: x**0 = 1 */
-	if((iy|ly)==0) return one; 	
+	if((iy|ly)==0) return one;
 
     /* +-NaN return x+y */
 	if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
-	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) 
-		return x+y;	
+	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
+		return x+y;
 
     /* determine if y is an odd int when x < 0
      * yisint = 0	... y is not an integer
@@ -130,7 +130,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
      * yisint = 2	... y is an even int
      */
 	yisint  = 0;
-	if(hx<0) {	
+	if(hx<0) {
 	    if(iy>=0x43400000) yisint = 2; /* even integer y */
 	    else if(iy>=0x3ff00000) {
 		k = (iy>>20)-0x3ff;	   /* exponent */
@@ -141,11 +141,11 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 		    j = iy>>(20-k);
 		    if((j<<(20-k))==iy) yisint = 2-(j&1);
 		}
-	    }		
-	} 
+	    }
+	}
 
     /* special value of y */
-	if(ly==0) { 	
+	if(ly==0) {
 	    if (iy==0x7ff00000) {	/* y is +-inf */
 	        if(((ix-0x3ff00000)|lx)==0)
 		    return  y - y;	/* inf**+-1 is NaN */
@@ -153,14 +153,14 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 		    return (hy>=0)? y: zero;
 	        else			/* (|x|<1)**-,+inf = inf,0 */
 		    return (hy<0)?-y: zero;
-	    } 
+	    }
 	    if(iy==0x3ff00000) {	/* y is  +-1 */
 		if(hy<0) return one/x; else return x;
 	    }
 	    if(hy==0x40000000) return x*x; /* y is  2 */
 	    if(hy==0x3fe00000) {	/* y is  0.5 */
 		if(hx>=0)	/* x >= +0 */
-		return __ieee754_sqrt(x);	
+		return __ieee754_sqrt(x);
 	    }
 	}
 
@@ -173,13 +173,13 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 		if(hx<0) {
 		    if(((ix-0x3ff00000)|yisint)==0) {
 			z = (z-z)/(z-z); /* (-1)**non-int is NaN */
-		    } else if(yisint==1) 
+		    } else if(yisint==1)
 			z = -z;		/* (x<0)**odd = -(|x|**odd) */
 		}
 		return z;
 	    }
 	}
-    
+
     /* (x<0)**(non-int) is NaN */
 	if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
 
@@ -192,7 +192,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	/* over/underflow if x is not close to one */
 	    if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
 	    if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
-	/* now |1-x| is tiny <= 2**-20, suffice to compute 
+	/* now |1-x| is tiny <= 2**-20, suffice to compute
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
 	    t = x-1;		/* t has 20 trailing zeros */
 	    w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
@@ -289,7 +289,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	    n = ((n&0x000fffff)|0x00100000)>>(20-k);
 	    if(j<0) n = -n;
 	    p_h -= t;
-	} 
+	}
 	t = p_l+p_h;
 	SET_LOW_WORD(t,0);
 	u = t*lg2_h;

libm/e_rem_pio2.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,8 +15,8 @@ static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $
 #endif
 
 /* __ieee754_rem_pio2(x,y)
- * 
- * return the remainder of x rem pi/2 in y[0]+y[1] 
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
  * use __kernel_rem_pio2()
  */
 
@@ -24,24 +24,24 @@ static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $
 #include "math_private.h"
 
 /*
- * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 
+ * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
  */
 #ifdef __STDC__
 static const int32_t two_over_pi[] = {
 #else
 static int32_t two_over_pi[] = {
 #endif
-0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
-0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
-0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
-0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
-0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
-0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
-0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
-0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
-0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
-0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
-0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
+0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
+0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
+0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
+0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
+0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
+0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
+0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
+0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
+0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
+0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
+0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
 };
 
 #ifdef __STDC__
@@ -68,9 +68,9 @@ static int32_t npio2_hw[] = {
  */
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
@@ -100,7 +100,7 @@ pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 	if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
 	    {y[0] = x; y[1] = 0; return 0;}
 	if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
-	    if(hx>0) { 
+	    if(hx>0) {
 		z = x - pio2_1;
 		if(ix!=0x3ff921fb) { 	/* 33+53 bit pi is good enough */
 		    y[0] = z - pio2_1t;
@@ -130,27 +130,27 @@ pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 	    fn = (double)n;
 	    r  = t-fn*pio2_1;
 	    w  = fn*pio2_1t;	/* 1st round good to 85 bit */
-	    if(n<32&&ix!=npio2_hw[n-1]) {	
+	    if(n<32&&ix!=npio2_hw[n-1]) {
 		y[0] = r-w;	/* quick check no cancellation */
 	    } else {
 	        u_int32_t high;
 	        j  = ix>>20;
-	        y[0] = r-w; 
+	        y[0] = r-w;
 		GET_HIGH_WORD(high,y[0]);
 	        i = j-((high>>20)&0x7ff);
 	        if(i>16) {  /* 2nd iteration needed, good to 118 */
 		    t  = r;
-		    w  = fn*pio2_2;	
+		    w  = fn*pio2_2;
 		    r  = t-w;
-		    w  = fn*pio2_2t-((t-r)-w);	
+		    w  = fn*pio2_2t-((t-r)-w);
 		    y[0] = r-w;
 		    GET_HIGH_WORD(high,y[0]);
 		    i = j-((high>>20)&0x7ff);
 		    if(i>49)  {	/* 3rd iteration need, 151 bits acc */
 		    	t  = r;	/* will cover all possible cases */
-		    	w  = fn*pio2_3;	
+		    	w  = fn*pio2_3;
 		    	r  = t-w;
-		    	w  = fn*pio2_3t-((t-r)-w);	
+		    	w  = fn*pio2_3t-((t-r)-w);
 		    	y[0] = r-w;
 		    }
 		}
@@ -159,7 +159,7 @@ pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
 	    if(hx<0) 	{y[0] = -y[0]; y[1] = -y[1]; return -n;}
 	    else	 return n;
 	}
-    /* 
+    /*
      * all other (large) arguments
      */
 	if(ix>=0x7ff00000) {		/* x is inf or NaN */

libm/e_remainder.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,11 +15,11 @@ static char rcsid[] = "$NetBSD: e_remainder.c,v 1.8 1995/05/10 20:46:05 jtc Exp
 #endif
 
 /* __ieee754_remainder(x,p)
- * Return :                  
- * 	returns  x REM p  =  x - [x/p]*p as if in infinite 
- * 	precise arithmetic, where [x/p] is the (infinite bit) 
+ * Return :
+ * 	returns  x REM p  =  x - [x/p]*p as if in infinite
+ * 	precise arithmetic, where [x/p] is the (infinite bit)
  *	integer nearest x/p (in half way case choose the even one).
- * Method : 
+ * Method :
  *	Based on fmod() return x-[x/p]chopped*p exactlp.
  */
 

libm/e_scalb.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -16,7 +16,7 @@ static char rcsid[] = "$NetBSD: e_scalb.c,v 1.6 1995/05/10 20:46:09 jtc Exp $";
 
 /*
  * __ieee754_scalb(x, fn) is provide for
- * passing various standard test suite. One 
+ * passing various standard test suite. One
  * should use scalbn() instead.
  */
 

libm/e_sinh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,15 +15,15 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
 #endif
 
 /* __ieee754_sinh(x)
- * Method : 
+ * Method :
  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
- *	1. Replace x by |x| (sinh(-x) = -sinh(x)). 
- *	2. 
+ *	1. Replace x by |x| (sinh(-x) = -sinh(x)).
+ *	2.
  *		                                    E + E/(E+1)
  *	    0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)
  *			       			        2
  *
- *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2 
+ *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2
  *	    lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)
  *	    ln2ovft  <  x	    :  sinh(x) := x*shuge (overflow)
  *
@@ -47,7 +47,7 @@ static double one = 1.0, shuge = 1.0e307;
 	double __ieee754_sinh(x)
 	double x;
 #endif
-{	
+{
 	double t,w,h;
 	int32_t ix,jx;
 	u_int32_t lx;
@@ -57,7 +57,7 @@ static double one = 1.0, shuge = 1.0e307;
 	ix = jx&0x7fffffff;
 
     /* x is INF or NaN */
-	if(ix>=0x7ff00000) return x+x;	
+	if(ix>=0x7ff00000) return x+x;
 
 	h = 0.5;
 	if (jx<0) h = -h;

libm/e_sqrt.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -19,10 +19,10 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
  *           ------------------------------------------
  *	     |  Use the hardware sqrt if you have one |
  *           ------------------------------------------
- * Method: 
- *   Bit by bit method using integer arithmetic. (Slow, but portable) 
+ * Method:
+ *   Bit by bit method using integer arithmetic. (Slow, but portable)
  *   1. Normalization
- *	Scale x to y in [1,4) with even powers of 2: 
+ *	Scale x to y in [1,4) with even powers of 2:
  *	find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
  *		sqrt(x) = 2^k * sqrt(y)
  *   2. Bit by bit computation
@@ -31,9 +31,9 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
  *                                     i+1         2
  *	    s  = 2*q , and	y  =  2   * ( y - q  ).		(1)
  *	     i      i            i                 i
- *                                                        
- *	To compute q    from q , one checks whether 
- *		    i+1       i                       
+ *
+ *	To compute q    from q , one checks whether
+ *		    i+1       i
  *
  *			      -(i+1) 2
  *			(q + 2      ) <= y.			(2)
@@ -43,12 +43,12 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
  *		 	       i+1   i             i+1   i
  *
  *	With some algebric manipulation, it is not difficult to see
- *	that (2) is equivalent to 
+ *	that (2) is equivalent to
  *                             -(i+1)
  *			s  +  2       <= y			(3)
  *			 i                i
  *
- *	The advantage of (3) is that s  and y  can be computed by 
+ *	The advantage of (3) is that s  and y  can be computed by
  *				      i      i
  *	the following recurrence formula:
  *	    if (3) is false
@@ -60,10 +60,10 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
  *                         -i                     -(i+1)
  *	    s	  =  s  + 2  ,  y    = y  -  s  - 2  		(5)
  *           i+1      i          i+1    i     i
- *				
- *	One may easily use induction to prove (4) and (5). 
+ *
+ *	One may easily use induction to prove (4) and (5).
  *	Note. Since the left hand side of (3) contain only i+2 bits,
- *	      it does not necessary to do a full (53-bit) comparison 
+ *	      it does not necessary to do a full (53-bit) comparison
  *	      in (3).
  *   3. Final rounding
  *	After generating the 53 bits result, we compute one more bit.
@@ -73,7 +73,7 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
  *	The rounding mode can be detected by checking whether
  *	huge + tiny is equal to huge, and whether huge - tiny is
  *	equal to huge for some floating point number "huge" and "tiny".
- *		
+ *
  * Special cases:
  *	sqrt(+-0) = +-0 	... exact
  *	sqrt(inf) = inf
@@ -101,17 +101,17 @@ static	double	one	= 1.0, tiny=1.0e-300;
 #endif
 {
 	double z;
-	int32_t sign = (int)0x80000000; 
+	int32_t sign = (int)0x80000000;
 	int32_t ix0,s0,q,m,t,i;
 	u_int32_t r,t1,s1,ix1,q1;
 
 	EXTRACT_WORDS(ix0,ix1,x);
 
     /* take care of Inf and NaN */
-	if((ix0&0x7ff00000)==0x7ff00000) {			
+	if((ix0&0x7ff00000)==0x7ff00000) {
 	    return x*x+x;		/* sqrt(NaN)=NaN, sqrt(+inf)=+inf
 					   sqrt(-inf)=sNaN */
-	} 
+	}
     /* take care of zero */
 	if(ix0<=0) {
 	    if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
@@ -145,12 +145,12 @@ static	double	one	= 1.0, tiny=1.0e-300;
 	r = 0x00200000;		/* r = moving bit from right to left */
 
 	while(r!=0) {
-	    t = s0+r; 
-	    if(t<=ix0) { 
-		s0   = t+r; 
-		ix0 -= t; 
-		q   += r; 
-	    } 
+	    t = s0+r;
+	    if(t<=ix0) {
+		s0   = t+r;
+		ix0 -= t;
+		q   += r;
+	    }
 	    ix0 += ix0 + ((ix1&sign)>>31);
 	    ix1 += ix1;
 	    r>>=1;
@@ -158,9 +158,9 @@ static	double	one	= 1.0, tiny=1.0e-300;
 
 	r = sign;
 	while(r!=0) {
-	    t1 = s1+r; 
+	    t1 = s1+r;
 	    t  = s0;
-	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 
+	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
 		s1  = t1+r;
 		if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
 		ix0 -= t;
@@ -181,7 +181,7 @@ static	double	one	= 1.0, tiny=1.0e-300;
 	        if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
 		else if (z>one) {
 		    if (q1==(u_int32_t)0xfffffffe) q+=1;
-		    q1+=2; 
+		    q1+=2;
 		} else
 	            q1 += (q1&1);
 	    }
@@ -197,18 +197,18 @@ static	double	one	= 1.0, tiny=1.0e-300;
 /*
 Other methods  (use floating-point arithmetic)
 -------------
-(This is a copy of a drafted paper by Prof W. Kahan 
+(This is a copy of a drafted paper by Prof W. Kahan
 and K.C. Ng, written in May, 1986)
 
-	Two algorithms are given here to implement sqrt(x) 
+	Two algorithms are given here to implement sqrt(x)
 	(IEEE double precision arithmetic) in software.
 	Both supply sqrt(x) correctly rounded. The first algorithm (in
 	Section A) uses newton iterations and involves four divisions.
 	The second one uses reciproot iterations to avoid division, but
 	requires more multiplications. Both algorithms need the ability
-	to chop results of arithmetic operations instead of round them, 
+	to chop results of arithmetic operations instead of round them,
 	and the INEXACT flag to indicate when an arithmetic operation
-	is executed exactly with no roundoff error, all part of the 
+	is executed exactly with no roundoff error, all part of the
 	standard (IEEE 754-1985). The ability to perform shift, add,
 	subtract and logical AND operations upon 32-bit words is needed
 	too, though not part of the standard.
@@ -218,7 +218,7 @@ A.  sqrt(x) by Newton Iteration
    (1)	Initial approximation
 
 	Let x0 and x1 be the leading and the trailing 32-bit words of
-	a floating point number x (in IEEE double format) respectively 
+	a floating point number x (in IEEE double format) respectively
 
 	    1    11		     52				  ...widths
 	   ------------------------------------------------------
@@ -226,7 +226,7 @@ A.  sqrt(x) by Newton Iteration
 	   ------------------------------------------------------
 	      msb    lsb  msb				      lsb ...order
 
- 
+
 	     ------------------------  	     ------------------------
 	x0:  |s|   e    |    f1     |	 x1: |          f2           |
 	     ------------------------  	     ------------------------
@@ -251,7 +251,7 @@ A.  sqrt(x) by Newton Iteration
 
     (2)	Iterative refinement
 
-	Apply Heron's rule three times to y, we have y approximates 
+	Apply Heron's rule three times to y, we have y approximates
 	sqrt(x) to within 1 ulp (Unit in the Last Place):
 
 		y := (y+x/y)/2		... almost 17 sig. bits
@@ -276,12 +276,12 @@ A.  sqrt(x) by Newton Iteration
 	it requires more multiplications and additions. Also x must be
 	scaled in advance to avoid spurious overflow in evaluating the
 	expression 3y*y+x. Hence it is not recommended uless division
-	is slow. If division is very slow, then one should use the 
+	is slow. If division is very slow, then one should use the
 	reciproot algorithm given in section B.
 
     (3) Final adjustment
 
-	By twiddling y's last bit it is possible to force y to be 
+	By twiddling y's last bit it is possible to force y to be
 	correctly rounded according to the prevailing rounding mode
 	as follows. Let r and i be copies of the rounding mode and
 	inexact flag before entering the square root program. Also we
@@ -312,7 +312,7 @@ A.  sqrt(x) by Newton Iteration
 	        I := i;	 		... restore inexact flag
 	        R := r;  		... restore rounded mode
 	        return sqrt(x):=y.
-		    
+
     (4)	Special cases
 
 	Square root of +inf, +-0, or NaN is itself;
@@ -331,7 +331,7 @@ B.  sqrt(x) by Reciproot Iteration
 	    k := 0x5fe80000 - (x0>>1);
 	    y0:= k - T2[63&(k>>14)].	... y ~ 1/sqrt(x) to 7.8 bits
 
-	Here k is a 32-bit integer and T2[] is an integer array 
+	Here k is a 32-bit integer and T2[] is an integer array
 	containing correction terms. Now magically the floating
 	value of y (y's leading 32-bit word is y0, the value of
 	its trailing word y1 is set to zero) approximates 1/sqrt(x)
@@ -352,9 +352,9 @@ B.  sqrt(x) by Reciproot Iteration
 
 	Apply Reciproot iteration three times to y and multiply the
 	result by x to get an approximation z that matches sqrt(x)
-	to about 1 ulp. To be exact, we will have 
+	to about 1 ulp. To be exact, we will have
 		-1ulp < sqrt(x)-z<1.0625ulp.
-	
+
 	... set rounding mode to Round-to-nearest
 	   y := y*(1.5-0.5*x*y*y)	... almost 15 sig. bits to 1/sqrt(x)
 	   y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
@@ -363,14 +363,14 @@ B.  sqrt(x) by Reciproot Iteration
 	   z := z + 0.5*z*(1-z*y)	... about 1 ulp to sqrt(x)
 
 	Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
-	(a) the term z*y in the final iteration is always less than 1; 
+	(a) the term z*y in the final iteration is always less than 1;
 	(b) the error in the final result is biased upward so that
 		-1 ulp < sqrt(x) - z < 1.0625 ulp
 	    instead of |sqrt(x)-z|<1.03125ulp.
 
     (3)	Final adjustment
 
-	By twiddling y's last bit it is possible to force y to be 
+	By twiddling y's last bit it is possible to force y to be
 	correctly rounded according to the prevailing rounding mode
 	as follows. Let r and i be copies of the rounding mode and
 	inexact flag before entering the square root program. Also we
@@ -410,27 +410,27 @@ B.  sqrt(x) by Reciproot Iteration
 	    I := 1;		... Raise Inexact flag: z is not exact
 	else {
 	    j := 1 - [(x0>>20)&1]	... j = logb(x) mod 2
-	    k := z1 >> 26;		... get z's 25-th and 26-th 
+	    k := z1 >> 26;		... get z's 25-th and 26-th
 					    fraction bits
 	    I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
 	}
 	R:= r		... restore rounded mode
 	return sqrt(x):=z.
 
-	If multiplication is cheaper then the foregoing red tape, the 
+	If multiplication is cheaper then the foregoing red tape, the
 	Inexact flag can be evaluated by
 
 	    I := i;
 	    I := (z*z!=x) or I.
 
-	Note that z*z can overwrite I; this value must be sensed if it is 
+	Note that z*z can overwrite I; this value must be sensed if it is
 	True.
 
 	Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
 	zero.
 
 		    --------------------
-		z1: |        f2        | 
+		z1: |        f2        |
 		    --------------------
 		bit 31		   bit 0
 
@@ -447,7 +447,7 @@ B.  sqrt(x) by Reciproot Iteration
 	11			01		even
 	-------------------------------------------------
 
-    (4)	Special cases (see (4) of Section A).	
- 
+    (4)	Special cases (see (4) of Section A).
+
  */
- 
+

libm/fp_private.h

@@ -56,13 +56,13 @@
 
 //#define    Infinity(x)       ( x.hex.exponent & ExpMask ) == ExpMask
 #define      dInfinity(x)      ( x.hex.high & dExpMask ) == dExpMask
-#define      sInfinity(x)      ( ( x.hexsgl << 1 ) & sExpMask ) == sExpMask      
+#define      sInfinity(x)      ( ( x.hexsgl << 1 ) & sExpMask ) == sExpMask
 
 //#define    Exponent(x)       x.hex.exponent & ExpMask
 #define      dExponent(x)      x.hex.high & dExpMask
 #define      sExponent(x)      ( ( x.hexsgl << 1 ) & sExpMask )
 
-#define      sZero(x)          ( x.hexsgl & sSgnMask ) == 0 
+#define      sZero(x)          ( x.hexsgl & sSgnMask ) == 0
 //#define    Sign(x)           ( x.hex.exponent & SgnMask ) == SgnMask
 
 /*******************************************************************************

libm/fpmacros.c

@@ -1,6 +1,6 @@
 /***********************************************************************
 **  File:  fpmacros.c
-**   
+**
 **  Contains:  C source code for implementations of floating-point
 **             functions which involve float format numbers, as
 **             defined in header <fp.h>.  In particular, this file
@@ -8,16 +8,16 @@
 **              __fpclassify(d,f), __isnormal(d,f), __isfinite(d,f),
 **             __isnan(d,f), and __signbit(d,f).  This file targets
 **             PowerPC platforms.
-**            
+**
 **  Written by:   Robert A. Murley, Ali Sazegari
-**   
+**
 **  Copyright:   c 2001 by Apple Computer, Inc., all rights reserved
-**   
+**
 **  Change History (most recent first):
 **
-**     07 Jul 01   ram      First created from fpfloatfunc.c, fp.c, 
+**     07 Jul 01   ram      First created from fpfloatfunc.c, fp.c,
 **				classify.c and sign.c in MathLib v3 Mac OS9.
-**            
+**
 ***********************************************************************/
 
 #include     <features.h>
@@ -33,48 +33,48 @@
 /***********************************************************************
    int __fpclassifyf(float x) returns the classification code of the
    argument x, as defined in <fp.h>.
-   
+
    Exceptions:  INVALID signaled if x is a signaling NaN; in this case,
                 the FP_QNAN code is returned.
-   
+
    Calls:  none
 ***********************************************************************/
 
 int __fpclassifyf ( float x )
 {
    unsigned int iexp;
-   
+
    union {
       u_int32_t lval;
       float fval;
    } z;
-   
+
    z.fval = x;
-   iexp = z.lval & FEXP_MASK;                 /* isolate float exponent */ 
-   
+   iexp = z.lval & FEXP_MASK;                 /* isolate float exponent */
+
    if (iexp == FEXP_MASK) {                   /* NaN or INF case */
       if ((z.lval & 0x007fffff) == 0)
          return FP_INFINITE;
 	return FP_NAN;
    }
-   
+
    if (iexp != 0)                             /* normal float */
       return FP_NORMAL;
-      
+
    if (x == 0.0)
       return FP_ZERO;             /* zero */
    else
       return FP_SUBNORMAL;        /* must be subnormal */
 }
-   
+
 
 /***********************************************************************
-      Function __fpclassify,                                                 
-      Implementation of classify of a double number for the PowerPC.          
-                                                                              
+      Function __fpclassify,
+      Implementation of classify of a double number for the PowerPC.
+
    Exceptions:  INVALID signaled if x is a signaling NaN; in this case,
                 the FP_QNAN code is returned.
-   
+
    Calls:  none
 ***********************************************************************/
 
@@ -86,16 +86,16 @@ int __fpclassify ( double arg )
             dHexParts hex;
             double dbl;
             } x;
-      
+
 	x.dbl = arg;
-	
+
 	exponent = x.hex.high & dExpMask;
 	if ( exponent == dExpMask )
 		{
 		if ( ( ( x.hex.high & dHighMan ) | x.hex.low ) == 0 )
 			return FP_INFINITE;
 		else
-            	return FP_NAN; 
+            	return FP_NAN;
 		}
 	else if ( exponent != 0)
 		return FP_NORMAL;
@@ -111,10 +111,10 @@ int __fpclassify ( double arg )
 /***********************************************************************
    int __isnormalf(float x) returns nonzero if and only if x is a
    normalized float number and zero otherwise.
-   
+
    Exceptions:  INVALID is raised if x is a signaling NaN; in this case,
                 zero is returned.
-   
+
    Calls:  none
 ***********************************************************************/
 
@@ -125,44 +125,44 @@ int __isnormalf ( float x )
       u_int32_t lval;
       float fval;
    } z;
-   
+
    z.fval = x;
    iexp = z.lval & FEXP_MASK;                 /* isolate float exponent */
    return ((iexp != FEXP_MASK) && (iexp != 0));
 }
-   
+
 
 int __isnormal ( double x )
 {
-	return ( __fpclassify ( x ) == FP_NORMAL ); 
+	return ( __fpclassify ( x ) == FP_NORMAL );
 }
 
 
 /***********************************************************************
    int __isfinitef(float x) returns nonzero if and only if x is a
    finite (normal, subnormal, or zero) float number and zero otherwise.
-   
+
    Exceptions:  INVALID is raised if x is a signaling NaN; in this case,
                 zero is returned.
-   
+
    Calls:  none
 ***********************************************************************/
 
 int __finitef ( float x )
-{   
+{
    union {
       u_int32_t lval;
       float fval;
    } z;
-   
+
    z.fval = x;
    return ((z.lval & FEXP_MASK) != FEXP_MASK);
 }
 weak_alias (__finitef, finitef)
-   
+
 int __finite ( double x )
 {
-	return ( __fpclassify ( x ) >= FP_ZERO ); 
+	return ( __fpclassify ( x ) >= FP_ZERO );
 }
 weak_alias (__finite, finite)
 
@@ -170,28 +170,28 @@ weak_alias (__finite, finite)
 /***********************************************************************
    int __signbitf(float x) returns nonzero if and only if the sign
    bit of x is set and zero otherwise.
-   
+
    Exceptions:  INVALID is raised if x is a signaling NaN.
-   
+
    Calls:  none
 ***********************************************************************/
 
 int __signbitf ( float x )
-{   
+{
    union {
       u_int32_t lval;
       float fval;
    } z;
-   
+
    z.fval = x;
    return ((z.lval & SIGN_MASK) != 0);
 }
 
 
 /***********************************************************************
-      Function sign of a double.                                              
-      Implementation of sign bit for the PowerPC.                             
-   
+      Function sign of a double.
+      Implementation of sign bit for the PowerPC.
+
    Calls:  none
 ***********************************************************************/
 
@@ -252,20 +252,20 @@ weak_alias (__isinfl, isinfl);
 /***********************************************************************
    int __isnanf(float x) returns nonzero if and only if x is a
    NaN and zero otherwise.
-   
+
    Exceptions:  INVALID is raised if x is a signaling NaN; in this case,
                 nonzero is returned.
-   
+
    Calls:  none
 ***********************************************************************/
 
 int __isnanf ( float x )
-{   
+{
    union {
       u_int32_t lval;
       float fval;
    } z;
-   
+
    z.fval = x;
    return (((z.lval&FEXP_MASK) == FEXP_MASK) && ((z.lval&FFRAC_MASK) != 0));
 }
@@ -274,7 +274,7 @@ weak_alias (__isnanf, isnanf);
 int __isnan ( double x )
 {
 	int class = __fpclassify(x);
-	return ( class == FP_NAN ); 
+	return ( class == FP_NAN );
 }
 weak_alias (__isnan, isnan);
 
@@ -282,7 +282,7 @@ weak_alias (__isnan, isnan);
 int __isnanl ( long double x )
 {
 	int class = __fpclassify(x);
-	return ( class == FP_NAN ); 
+	return ( class == FP_NAN );
 }
 weak_alias (__isnanl, isnanl);
 #endif

libm/k_cos.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -18,7 +18,7 @@ static char rcsid[] = "$NetBSD: k_cos.c,v 1.8 1995/05/10 20:46:22 jtc Exp $";
  * __kernel_cos( x,  y )
  * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
  * Input x is assumed to be bounded by ~pi/4 in magnitude.
- * Input y is the tail of x. 
+ * Input y is the tail of x.
  *
  * Algorithm
  *	1. Since cos(-x) = cos(x), we need only to consider positive x.
@@ -28,15 +28,15 @@ static char rcsid[] = "$NetBSD: k_cos.c,v 1.8 1995/05/10 20:46:22 jtc Exp $";
  *		  	                 4            14
  *	   	cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
  *	   where the remez error is
- *	
+ *
  * 	|              2     4     6     8     10    12     14 |     -58
  * 	|cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
- * 	|    					               | 
- * 
- * 	               4     6     8     10    12     14 
+ * 	|    					               |
+ *
+ * 	               4     6     8     10    12     14
  *	4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
  *	       cos(x) = 1 - x*x/2 + r
- *	   since cos(x+y) ~ cos(x) - sin(x)*y 
+ *	   since cos(x+y) ~ cos(x) - sin(x)*y
  *			  ~ cos(x) - x*y,
  *	   a correction term is necessary in cos(x) and hence
  *		cos(x+y) = 1 - (x*x/2 - (r - x*y))
@@ -53,9 +53,9 @@ static char rcsid[] = "$NetBSD: k_cos.c,v 1.8 1995/05/10 20:46:22 jtc Exp $";
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
@@ -81,7 +81,7 @@ C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
 	}
 	z  = x*x;
 	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
-	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */ 
+	if(ix < 0x3FD33333) 			/* if |x| < 0.3 */
 	    return one - (0.5*z - (z*r - x*y));
 	else {
 	    if(ix > 0x3fe90000) {		/* x > 0.78125 */

libm/k_rem_pio2.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -17,12 +17,12 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
 /*
  * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
  * double x[],y[]; int e0,nx,prec; int ipio2[];
- * 
- * __kernel_rem_pio2 return the last three digits of N with 
+ *
+ * __kernel_rem_pio2 return the last three digits of N with
  *		y = x - N*pi/2
  * so that |y| < pi/2.
  *
- * The method is to compute the integer (mod 8) and fraction parts of 
+ * The method is to compute the integer (mod 8) and fraction parts of
  * (2/pi)*x without doing the full multiplication. In general we
  * skip the part of the product that are known to be a huge integer (
  * more accurately, = 0 mod 8 ). Thus the number of operations are
@@ -31,10 +31,10 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
  * (2/pi) is represented by an array of 24-bit integers in ipio2[].
  *
  * Input parameters:
- * 	x[]	The input value (must be positive) is broken into nx 
+ * 	x[]	The input value (must be positive) is broken into nx
  *		pieces of 24-bit integers in double precision format.
- *		x[i] will be the i-th 24 bit of x. The scaled exponent 
- *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 
+ *		x[i] will be the i-th 24 bit of x. The scaled exponent
+ *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
  *		match x's up to 24 bits.
  *
  *		Example of breaking a double positive z into x[0]+x[1]+x[2]:
@@ -71,8 +71,8 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
  *			3	113 bits (quad)
  *
  *	ipio2[]
- *		integer array, contains the (24*i)-th to (24*i+23)-th 
- *		bit of 2/pi after binary point. The corresponding 
+ *		integer array, contains the (24*i)-th to (24*i+23)-th
+ *		bit of 2/pi after binary point. The corresponding
  *		floating value is
  *
  *			ipio2[i] * 2^(-24(i+1)).
@@ -87,8 +87,8 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
  *		in the computation. The recommended value is 2,3,4,
  *		6 for single, double, extended,and quad.
  *
- * 	jz	local integer variable indicating the number of 
- *		terms of ipio2[] used. 
+ * 	jz	local integer variable indicating the number of
+ *		terms of ipio2[] used.
  *
  *	jx	nx - 1
  *
@@ -108,9 +108,9 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
  *		exponent for q[i] would be q0-24*i.
  *
  *	PIo2[]	double precision array, obtained by cutting pi/2
- *		into 24 bits chunks. 
+ *		into 24 bits chunks.
  *
- *	f[]	ipio2[] in floating point 
+ *	f[]	ipio2[] in floating point
  *
  *	iq[]	integer array by breaking up q[] in 24-bits chunk.
  *
@@ -124,9 +124,9 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
 
 /*
  * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
- * compiler will convert from decimal to binary accurately enough 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
 
@@ -136,7 +136,7 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
 #ifdef __STDC__
 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
 #else
-static int init_jk[] = {2,3,4,6}; 
+static int init_jk[] = {2,3,4,6};
 #endif
 
 #ifdef __STDC__
@@ -155,9 +155,9 @@ static double PIo2[] = {
 };
 
 #ifdef __STDC__
-static const double			
+static const double
 #else
-static double			
+static double
 #endif
 zero   = 0.0,
 one    = 1.0,
@@ -165,9 +165,9 @@ two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
 
 #ifdef __STDC__
-	int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2) 
+	int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
 #else
-	int __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) 	
+	int __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
 	double x[], y[]; int e0,nx,prec; int32_t ipio2[];
 #endif
 {
@@ -211,7 +211,7 @@ recompute:
 	    i  = (iq[jz-1]>>(24-q0)); n += i;
 	    iq[jz-1] -= i<<(24-q0);
 	    ih = iq[jz-1]>>(23-q0);
-	} 
+	}
 	else if(q0==0) ih = iq[jz-1]>>23;
 	else if(z>=0.5) ih=2;
 
@@ -262,7 +262,7 @@ recompute:
 	    while(iq[jz]==0) { jz--; q0-=24;}
 	} else { /* break z into 24-bit if necessary */
 	    z = scalbn(z,-q0);
-	    if(z>=two24) { 
+	    if(z>=two24) {
 		fw = (double)((int32_t)(twon24*z));
 		iq[jz] = (int32_t)(z-two24*fw);
 		jz += 1; q0 += 24;
@@ -287,29 +287,29 @@ recompute:
 	    case 0:
 		fw = 0.0;
 		for (i=jz;i>=0;i--) fw += fq[i];
-		y[0] = (ih==0)? fw: -fw; 
+		y[0] = (ih==0)? fw: -fw;
 		break;
 	    case 1:
 	    case 2:
 		fw = 0.0;
-		for (i=jz;i>=0;i--) fw += fq[i]; 
-		y[0] = (ih==0)? fw: -fw; 
+		for (i=jz;i>=0;i--) fw += fq[i];
+		y[0] = (ih==0)? fw: -fw;
 		fw = fq[0]-fw;
 		for (i=1;i<=jz;i++) fw += fq[i];
-		y[1] = (ih==0)? fw: -fw; 
+		y[1] = (ih==0)? fw: -fw;
 		break;
 	    case 3:	/* painful */
 		for (i=jz;i>0;i--) {
-		    fw      = fq[i-1]+fq[i]; 
+		    fw      = fq[i-1]+fq[i];
 		    fq[i]  += fq[i-1]-fw;
 		    fq[i-1] = fw;
 		}
 		for (i=jz;i>1;i--) {
-		    fw      = fq[i-1]+fq[i]; 
+		    fw      = fq[i-1]+fq[i];
 		    fq[i]  += fq[i-1]-fw;
 		    fq[i-1] = fw;
 		}
-		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; 
+		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
 		if(ih==0) {
 		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
 		} else {

libm/k_sin.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -18,24 +18,24 @@ static char rcsid[] = "$NetBSD: k_sin.c,v 1.8 1995/05/10 20:46:31 jtc Exp $";
  * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
  * Input x is assumed to be bounded by ~pi/4 in magnitude.
  * Input y is the tail of x.
- * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). 
+ * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
  *
  * Algorithm
- *	1. Since sin(-x) = -sin(x), we need only to consider positive x. 
+ *	1. Since sin(-x) = -sin(x), we need only to consider positive x.
  *	2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
  *	3. sin(x) is approximated by a polynomial of degree 13 on
  *	   [0,pi/4]
  *		  	         3            13
  *	   	sin(x) ~ x + S1*x + ... + S6*x
  *	   where
- *	
+ *
  * 	|sin(x)         2     4     6     8     10     12  |     -58
  * 	|----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
- * 	|  x 					           | 
- * 
+ * 	|  x 					           |
+ *
  *	4. sin(x+y) = sin(x) + sin'(x')*y
  *		    ~ sin(x) + (1-x*x/2)*y
- *	   For better accuracy, let 
+ *	   For better accuracy, let
  *		     3      2      2      2      2
  *		r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
  *	   then                   3    2
@@ -46,9 +46,9 @@ static char rcsid[] = "$NetBSD: k_sin.c,v 1.8 1995/05/10 20:46:31 jtc Exp $";
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */

libm/k_standard.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -33,7 +33,7 @@ static const double zero = 0.0;	/* used as const */
 static double zero = 0.0;	/* used as const */
 #endif
 
-/* 
+/*
  * Standard conformance (non-IEEE) on exception cases.
  * Mapping:
  *	1 -- acos(|x|>1)
@@ -58,7 +58,7 @@ static double zero = 0.0;	/* used as const */
  *	20-- pow(0.0,0.0)
  *	21-- pow(x,y) overflow
  *	22-- pow(x,y) underflow
- *	23-- pow(0,negative) 
+ *	23-- pow(0,negative)
  *	24-- pow(neg,non-integral)
  *	25-- sinh(finite) overflow
  *	26-- sqrt(negative)
@@ -82,14 +82,14 @@ static double zero = 0.0;	/* used as const */
 
 
 #ifdef __STDC__
-	double __kernel_standard(double x, double y, int type) 
+	double __kernel_standard(double x, double y, int type)
 #else
-	double __kernel_standard(x,y,type) 
+	double __kernel_standard(x,y,type)
 	double x,y; int type;
 #endif
 {
 	struct exception exc;
-#ifndef HUGE_VAL	/* this is the only routine that uses HUGE_VAL */ 
+#ifndef HUGE_VAL	/* this is the only routine that uses HUGE_VAL */
 #define HUGE_VAL inf
 	double inf = 0.0;
 
@@ -469,7 +469,7 @@ static double zero = 0.0;	/* used as const */
 		/* 0**neg */
 		exc.type = DOMAIN;
 		exc.name = type < 100 ? "pow" : "powf";
-		if (_LIB_VERSION == _SVID_) 
+		if (_LIB_VERSION == _SVID_)
 		  exc.retval = zero;
 		else
 		  exc.retval = -HUGE_VAL;
@@ -487,11 +487,11 @@ static double zero = 0.0;	/* used as const */
 		/* neg**non-integral */
 		exc.type = DOMAIN;
 		exc.name = type < 100 ? "pow" : "powf";
-		if (_LIB_VERSION == _SVID_) 
+		if (_LIB_VERSION == _SVID_)
 		    exc.retval = zero;
-		else 
+		else
 		    exc.retval = zero/zero;	/* X/Open allow NaN */
-		if (_LIB_VERSION == _POSIX_) 
+		if (_LIB_VERSION == _POSIX_)
 		   errno = EDOM;
 		else if (!matherr(&exc)) {
 		  if (_LIB_VERSION == _SVID_) {
@@ -649,7 +649,7 @@ static double zero = 0.0;	/* used as const */
                                 (void) WRITE2(": TLOSS error\n", 14);
                         }
                         errno = ERANGE;
-                }        
+                }
 		break;
 	    case 35:
 	    case 135:
@@ -665,7 +665,7 @@ static double zero = 0.0;	/* used as const */
                                 (void) WRITE2(": TLOSS error\n", 14);
                         }
                         errno = ERANGE;
-                }        
+                }
 		break;
 	    case 36:
 	    case 136:
@@ -681,7 +681,7 @@ static double zero = 0.0;	/* used as const */
                                 (void) WRITE2(": TLOSS error\n", 14);
                         }
                         errno = ERANGE;
-                }        
+                }
 		break;
 	    case 37:
 	    case 137:
@@ -697,7 +697,7 @@ static double zero = 0.0;	/* used as const */
                                 (void) WRITE2(": TLOSS error\n", 14);
                         }
                         errno = ERANGE;
-                }        
+                }
 		break;
 	    case 38:
 	    case 138:
@@ -713,7 +713,7 @@ static double zero = 0.0;	/* used as const */
                                 (void) WRITE2(": TLOSS error\n", 14);
                         }
                         errno = ERANGE;
-                }        
+                }
 		break;
 	    case 39:
 	    case 139:
@@ -729,7 +729,7 @@ static double zero = 0.0;	/* used as const */
                                 (void) WRITE2(": TLOSS error\n", 14);
                         }
                         errno = ERANGE;
-                }        
+                }
 		break;
 	    case 40:
 	    case 140:
@@ -778,5 +778,5 @@ static double zero = 0.0;	/* used as const */
 		}
 		break;
 	}
-	return exc.retval; 
+	return exc.retval;
 }

libm/k_tan.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -18,25 +18,25 @@ static char rcsid[] = "$NetBSD: k_tan.c,v 1.8 1995/05/10 20:46:37 jtc Exp $";
  * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
  * Input x is assumed to be bounded by ~pi/4 in magnitude.
  * Input y is the tail of x.
- * Input k indicates whether tan (if k=1) or 
+ * Input k indicates whether tan (if k=1) or
  * -1/tan (if k= -1) is returned.
  *
  * Algorithm
- *	1. Since tan(-x) = -tan(x), we need only to consider positive x. 
+ *	1. Since tan(-x) = -tan(x), we need only to consider positive x.
  *	2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
  *	3. tan(x) is approximated by a odd polynomial of degree 27 on
  *	   [0,0.67434]
  *		  	         3             27
  *	   	tan(x) ~ x + T1*x + ... + T13*x
  *	   where
- *	
+ *
  * 	        |tan(x)         2     4            26   |     -59.2
  * 	        |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
- * 	        |  x 					| 
- * 
+ * 	        |  x 					|
+ *
  *	   Note: tan(x+y) = tan(x) + tan'(x)*y
  *		          ~ tan(x) + (1+x*x)*y
- *	   Therefore, for better accuracy in computing tan(x+y), let 
+ *	   Therefore, for better accuracy in computing tan(x+y), let
  *		     3      2      2       2       2
  *		r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
  *	   then
@@ -51,9 +51,9 @@ static char rcsid[] = "$NetBSD: k_tan.c,v 1.8 1995/05/10 20:46:37 jtc Exp $";
 #include "math.h"
 #include "math_private.h"
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 one   =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 pio4  =  7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
@@ -116,7 +116,7 @@ T[] =  {
 	    return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r)));
 	}
 	if(iy==1) return w;
-	else {		/* if allow error up to 2 ulp, 
+	else {		/* if allow error up to 2 ulp,
 			   simply return -1.0/(x+r) here */
      /*  compute -1.0/(x+r) accurately */
 	    double a,t;

libm/math_private.h

@@ -4,7 +4,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -42,10 +42,10 @@
 
 #if (__BYTE_ORDER == __BIG_ENDIAN) || defined(__arm__) && !defined(__VFP_FP__)
 
-typedef union 
+typedef union
 {
   double value;
-  struct 
+  struct
   {
     u_int32_t msw;
     u_int32_t lsw;
@@ -56,10 +56,10 @@ typedef union
 
 #if (__BYTE_ORDER == __LITTLE_ENDIAN) && (!defined(__arm__) || defined(__VFP_FP__))
 
-typedef union 
+typedef union
 {
   double value;
-  struct 
+  struct
   {
     u_int32_t lsw;
     u_int32_t msw;
@@ -154,13 +154,13 @@ do {								\
 } while (0)
 
 /* ieee style elementary functions */
-extern double __ieee754_sqrt __P((double));			
-extern double __ieee754_acos __P((double));			
-extern double __ieee754_acosh __P((double));			
-extern double __ieee754_log __P((double));			
-extern double __ieee754_atanh __P((double));			
-extern double __ieee754_asin __P((double));			
-extern double __ieee754_atan2 __P((double,double));			
+extern double __ieee754_sqrt __P((double));
+extern double __ieee754_acos __P((double));
+extern double __ieee754_acosh __P((double));
+extern double __ieee754_log __P((double));
+extern double __ieee754_atanh __P((double));
+extern double __ieee754_asin __P((double));
+extern double __ieee754_atan2 __P((double,double));
 extern double __ieee754_exp __P((double));
 extern double __ieee754_cosh __P((double));
 extern double __ieee754_fmod __P((double,double));
@@ -187,7 +187,7 @@ extern double __ieee754_scalb __P((double,double));
 #endif
 
 /* fdlibm kernel function */
-extern double __kernel_standard __P((double,double,int));	
+extern double __kernel_standard __P((double,double,int));
 extern double __kernel_sin __P((double,double,int));
 extern double __kernel_cos __P((double,double));
 extern double __kernel_tan __P((double,double,int));
@@ -195,13 +195,13 @@ extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const int*));
 
 
 /* ieee style elementary float functions */
-extern float __ieee754_sqrtf __P((float));			
-extern float __ieee754_acosf __P((float));			
-extern float __ieee754_acoshf __P((float));			
-extern float __ieee754_logf __P((float));			
-extern float __ieee754_atanhf __P((float));			
-extern float __ieee754_asinf __P((float));			
-extern float __ieee754_atan2f __P((float,float));			
+extern float __ieee754_sqrtf __P((float));
+extern float __ieee754_acosf __P((float));
+extern float __ieee754_acoshf __P((float));
+extern float __ieee754_logf __P((float));
+extern float __ieee754_atanhf __P((float));
+extern float __ieee754_asinf __P((float));
+extern float __ieee754_atan2f __P((float,float));
 extern float __ieee754_expf __P((float));
 extern float __ieee754_coshf __P((float));
 extern float __ieee754_fmodf __P((float,float));

libm/powerpc/Makefile

@@ -6,14 +6,14 @@
 # math library for Apple's MacOS X/Darwin math library, which was
 # itself swiped from FreeBSD.  The original copyright information
 # is as follows:
-# 
+#
 #     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.
-# 
+#
 # It has been ported to work with uClibc and generally behave
 # by Erik Andersen <andersen@codepoet.org>
 #
@@ -66,6 +66,6 @@ $(OBJ): Makefile
 tags:
 	ctags -R
 
-clean: 
+clean:
 	$(RM) *.[oa] *~ core
 

libm/powerpc/s_ceil.c

@@ -52,24 +52,24 @@ double ceil ( double x )
 	register double y;
 	register unsigned long int xhi;
 	register int target;
-	
+
 	xInHex.dbl = x;
 	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
 	target = ( xInHex.words.hi < signMask );
-	
-	if ( xhi < 0x43300000ul ) 
+
+	if ( xhi < 0x43300000ul )
 /*******************************************************************************
 *      Is |x| < 2.0^52?                                                        *
 *******************************************************************************/
 		{
-		if ( xhi < 0x3ff00000ul ) 
+		if ( xhi < 0x3ff00000ul )
 /*******************************************************************************
 *      Is |x| < 1.0?                                                           *
 *******************************************************************************/
 			{
 			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
 				return ( x );
-			else 
+			else
 				{			                // inexact case
 				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
 				OldEnvironment.words.lo |= 0x02000000ul;
@@ -83,7 +83,7 @@ double ceil ( double x )
 /*******************************************************************************
 *      Is 1.0 < |x| < 2.0^52?                                                  *
 *******************************************************************************/
-		if ( target ) 
+		if ( target )
 			{
 			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
 			if ( y < x )
@@ -91,8 +91,8 @@ double ceil ( double x )
 			else
 				return ( y );
 			}
-		
-		else 
+
+		else
 			{
 			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
 			if ( y < x )

libm/powerpc/s_copysign.c

@@ -45,12 +45,12 @@ double copysign ( double arg2, double arg1 )
 /*******************************************************************************
 *     No need to flush NaNs out.                                               *
 *******************************************************************************/
-      
+
       x.dbl = arg1;
       y.dbl = arg2;
-      
+
       y.hex.high = y.hex.high & 0x7FFFFFFF;
       y.hex.high = ( y.hex.high | ( x.hex.high & dSgnMask ) );
-      
+
       return y.dbl;
       }

libm/powerpc/s_floor.c

@@ -52,24 +52,24 @@ double floor ( double x )
 	register double y;
 	register unsigned long int xhi;
 	register long int target;
-	
+
 	xInHex.dbl = x;
 	xhi = xInHex.words.hi & 0x7fffffffUL;	  // xhi is the high half of |x|
 	target = ( xInHex.words.hi < signMask );
-	
-	if ( xhi < 0x43300000ul ) 
+
+	if ( xhi < 0x43300000ul )
 /*******************************************************************************
 *      Is |x| < 2.0^52?                                                        *
 *******************************************************************************/
 		{
-		if ( xhi < 0x3ff00000ul ) 
+		if ( xhi < 0x3ff00000ul )
 /*******************************************************************************
 *      Is |x| < 1.0?                                                           *
 *******************************************************************************/
 			{
 			if ( ( xhi | xInHex.words.lo ) == 0ul )  // zero x is exact case
 				return ( x );
-			else 
+			else
 				{			                // inexact case
 				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
 				OldEnvironment.words.lo |= 0x02000000ul;
@@ -83,7 +83,7 @@ double floor ( double x )
 /*******************************************************************************
 *      Is 1.0 < |x| < 2.0^52?                                                  *
 *******************************************************************************/
-		if ( target ) 
+		if ( target )
 			{
 			y = ( x + twoTo52 ) - twoTo52;          // round at binary pt.
 			if ( y > x )
@@ -91,8 +91,8 @@ double floor ( double x )
 			else
 				return ( y );
 			}
-		
-		else 
+
+		else
 			{
 			y = ( x - twoTo52 ) + twoTo52;          // round at binary pt.
 			if ( y > x )

libm/powerpc/s_frexp.c

@@ -31,12 +31,12 @@ typedef union
         unsigned long int hi;
         unsigned long int lo;
 #else
-        unsigned long int lo; 
-        unsigned long int hi; 
+        unsigned long int lo;
+        unsigned long int hi;
 #endif
       } words;
       double dbl;
-      } DblInHex;       
+      } DblInHex;
 
 double frexp ( double value, int *eptr )
       {
@@ -49,7 +49,7 @@ double frexp ( double value, int *eptr )
       *eptr = 0;
 	if ( valueHead >= 0x7ff00000 || ( valueHead | argument.words.lo ) == 0 )
 		return value;		// 0, inf, or NaN
-	
+
 	if ( valueHead < 0x00100000 )
 		{	// denorm
 		argument.dbl = two54 * value;

libm/powerpc/s_ldexp.c

@@ -29,18 +29,18 @@ typedef union
         unsigned long int hi;
         unsigned long int lo;
 #else
-        unsigned long int lo; 
-        unsigned long int hi; 
+        unsigned long int lo;
+        unsigned long int hi;
 #endif
       } words;
       double dbl;
-      } DblInHex;       
+      } DblInHex;
 
-double ldexp ( double value, int exp ) 
+double ldexp ( double value, int exp )
       {
-      if ( exp > SHRT_MAX ) 
+      if ( exp > SHRT_MAX )
             exp = SHRT_MAX;
-      else if ( exp < -SHRT_MAX ) 
+      else if ( exp < -SHRT_MAX )
             exp = -SHRT_MAX;
       return scalb ( value, exp  );
       }

libm/powerpc/s_logb.c

@@ -32,8 +32,8 @@
 *     Standard 754.                                                            *
 *******************************************************************************/
 
-typedef union           
-      { 
+typedef union
+      {
       struct {
 #if defined(__BIG_ENDIAN__)
         unsigned long int hi;
@@ -62,38 +62,38 @@ double logb (  double x  )
       {
       DblInHex xInHex;
       long int shiftedExp;
-      
+
       xInHex.dbl = x;
       shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
-      
-      if ( shiftedExp == 2047 ) 
+
+      if ( shiftedExp == 2047 )
             {                                            // NaN or INF
             if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) )
                   return x;                              // NaN or +INF return x
             else
                   return -x;                             // -INF returns +INF
             }
-      
+
       if ( shiftedExp != 0 )                             // normal number
             shiftedExp -= 1023;                          // unbias exponent
-      
-      else if ( x == 0.0 ) 
+
+      else if ( x == 0.0 )
             {                                            // zero
             xInHex.words.hi = 0x0UL;                      // return -infinity
             return (  minusInf.dbl  );
             }
-      
-      else 
+
+      else
             {                                            // subnormal number
             xInHex.dbl *= twoTo52;                       // scale up
             shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
             shiftedExp -= 1075;                          // unbias exponent
             }
-      
+
       if ( shiftedExp == 0 )                             // zero result
             return ( 0.0 );
-      
-      else 
+
+      else
             {                                            // nonzero result
             xInHex.dbl = klTod;
             xInHex.words.lo += shiftedExp;

libm/powerpc/s_modf.c

@@ -1,18 +1,18 @@
 /*******************************************************************************
 **      File:   rndint.c
-**      
+**
 **      Contains: C source code for implementations of floating-point
 **                functions which round to integral value or format, as
 **                defined in header <fp.h>.  In particular, this file
 **                contains implementations of functions rint, nearbyint,
 **                rinttol, round, roundtol, trunc, modf and modfl.  This file
 **                targets PowerPC or Power platforms.
-**                        
+**
 **      Written by: A. Sazegari, Apple AltiVec Group
 **	    Created originally by Jon Okada, Apple Numerics Group
-**      
+**
 **      Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
-**      
+**
 **      Change History (most recent first):
 **
 **      13 Jul 01  ram  replaced --setflm calls with inline assembly
@@ -21,7 +21,7 @@
 **				1. removed double_t, put in double for now.
 **				2. removed iclass from nearbyint.
 **				3. removed wrong comments intrunc.
-**				4. 
+**				4.
 **      13 May 97  ali  made performance improvements in rint, rinttol, roundtol
 **                      and trunc by folding some of the taligent ideas into this
 **                      implementation.  nearbyint is faster than the one in taligent,
@@ -33,7 +33,7 @@
 **      18 Feb 93  ali  Changed the return value of fenv functions
 **                      feclearexcept and feraiseexcept to their new
 **                      NCEG X3J11.1/93-001 definitions.
-**      16 Dec 92  JPO  Removed __itrunc implementation to a 
+**      16 Dec 92  JPO  Removed __itrunc implementation to a
 **                      separate file.
 **      15 Dec 92  JPO  Added __itrunc implementation and modified
 **                      rinttol to include conversion from double
@@ -41,7 +41,7 @@
 **                      call __itrunc.
 **      10 Dec 92  JPO  Added modf (double) implementation.
 **      04 Dec 92  JPO  First created.
-**                        
+**
 *******************************************************************************/
 
 #include <limits.h>
@@ -81,14 +81,14 @@ static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
 *     This function calls fabs and copysign.		                         *
 *                                                                              *
 *******************************************************************************/
-   
+
 double nearbyint ( double x )
       {
 	double y;
 	double OldEnvironment;
-      
+
 	y = twoTo52;
-	
+
 	asm ("mffs %0" : "=f" (OldEnvironment));	/* get the environement */
 
       if ( fabs ( x ) >= y )                          /* huge case is exact */
@@ -99,9 +99,9 @@ double nearbyint ( double x )
             y = copysign ( y, x );
 //	restore old flags
 	asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
-      return ( y );      
+      return ( y );
 	}
-      
+
 /*******************************************************************************
 *                                                                              *
 *     The function rinttol converts its double argument to integral value      *
@@ -120,38 +120,38 @@ long int rinttol ( double x )
       DblInHex argument, OldEnvironment;
       unsigned long int xHead;
       register long int target;
-      
+
       argument.dbl = x;
       target = ( argument.words.hi < signMask );        // flag positive sign
       xHead = argument.words.hi & 0x7ffffffful;         // high 32 bits of x
-      
-      if ( target ) 
+
+      if ( target )
 /*******************************************************************************
 *    Sign of x is positive.                                                    *
 *******************************************************************************/
             {
-            if ( xHead < 0x41dffffful ) 
+            if ( xHead < 0x41dffffful )
                   {                                    // x is safely in long range
                   y = ( x + twoTo52 ) - twoTo52;       // round at binary point
                   argument.dbl = y + doubleToLong;     // force result into argument.words.lo
                   return ( ( long ) argument.words.lo );
                   }
-            
+
 		asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
-            
-            if ( xHead > 0x41dffffful ) 
+
+            if ( xHead > 0x41dffffful )
                   {                                    // x is safely out of long range
                   OldEnvironment.words.lo |= SET_INVALID;
 			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
                   return ( LONG_MAX );
                   }
-            
+
 /*******************************************************************************
 *     x > 0.0 and may or may not be out of range of long.                      *
 *******************************************************************************/
 
             y = ( x + twoTo52 ) - twoTo52;             // do rounding
-            if ( y > ( double ) LONG_MAX ) 
+            if ( y > ( double ) LONG_MAX )
                   {                                    // out of range of long
                   OldEnvironment.words.lo |= SET_INVALID;
 			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
@@ -160,32 +160,32 @@ long int rinttol ( double x )
             argument.dbl = y + doubleToLong;           // in range
             return ( ( long ) argument.words.lo );      // return result & flags
             }
-      
+
 /*******************************************************************************
 *    Sign of x is negative.                                                    *
 *******************************************************************************/
-      if ( xHead < 0x41e00000ul ) 
+      if ( xHead < 0x41e00000ul )
             {                                          // x is safely in long range
             y = ( x - twoTo52 ) + twoTo52;
             argument.dbl = y + doubleToLong;
             return ( ( long ) argument.words.lo );
             }
-      
+
 	asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
-      
-      if ( xHead > 0x41e00000ul ) 
+
+      if ( xHead > 0x41e00000ul )
             {                                          // x is safely out of long range
             OldEnvironment.words.lo |= SET_INVALID;
 		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
             return ( LONG_MIN );
             }
-      
+
 /*******************************************************************************
 *    x < 0.0 and may or may not be out of range of long.                       *
 *******************************************************************************/
 
       y = ( x - twoTo52 ) + twoTo52;                   // do rounding
-      if ( y < ( double ) LONG_MIN ) 
+      if ( y < ( double ) LONG_MIN )
             {                                          // out of range of long
             OldEnvironment.words.lo |= SET_INVALID;
 		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
@@ -204,30 +204,30 @@ long int rinttol ( double x )
 *     return value is not equal to the operand.                                *
 *                                                                              *
 *******************************************************************************/
-   
+
 double round ( double x )
-      {      
+      {
       DblInHex argument, OldEnvironment;
       register double y, z;
       register unsigned long int xHead;
       register long int target;
-      
+
       argument.dbl = x;
       xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
       target = ( argument.words.hi < signMask );     // flag positive sign
-      
-      if ( xHead < 0x43300000ul ) 
+
+      if ( xHead < 0x43300000ul )
 /*******************************************************************************
 *     Is |x| < 2.0^52?                                                        *
 *******************************************************************************/
             {
-            if ( xHead < 0x3ff00000ul ) 
+            if ( xHead < 0x3ff00000ul )
 /*******************************************************************************
 *     Is |x| < 1.0?                                                           *
 *******************************************************************************/
                   {
 			asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
-                  if ( xHead < 0x3fe00000ul ) 
+                  if ( xHead < 0x3fe00000ul )
 /*******************************************************************************
 *     Is |x| < 0.5?                                                           *
 *******************************************************************************/
@@ -235,7 +235,7 @@ double round ( double x )
                         if ( ( xHead | argument.words.lo ) != 0ul )
                               OldEnvironment.words.lo |= 0x02000000ul;
 				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
-                        if ( target ) 
+                        if ( target )
                               return ( 0.0 );
                         else
                               return ( -0.0 );
@@ -253,7 +253,7 @@ double round ( double x )
 /*******************************************************************************
 *     Is 1.0 < |x| < 2.0^52?                                                   *
 *******************************************************************************/
-            if ( target ) 
+            if ( target )
                   {                                     // positive x
                   y = ( x + twoTo52 ) - twoTo52;        // round at binary point
                   if ( y == x )                         // exact case
@@ -265,11 +265,11 @@ double round ( double x )
                   else
                         return ( y );
                   }
-            
+
 /*******************************************************************************
 *     Is x < 0?                                                                *
 *******************************************************************************/
-            else 
+            else
                   {
                   y = ( x - twoTo52 ) + twoTo52;        // round at binary point
                   if ( y == x )
@@ -302,19 +302,19 @@ double round ( double x )
 *******************************************************************************/
 
 long int roundtol ( double x )
-	{	
+	{
 	register double y, z;
 	DblInHex argument, OldEnvironment;
 	register unsigned long int xhi;
 	register long int target;
 	const DblInHex kTZ = {{ 0x0, 0x1 }};
 	const DblInHex kUP = {{ 0x0, 0x2 }};
-	
+
 	argument.dbl = x;
 	xhi = argument.words.hi & 0x7ffffffful;	        	// high 32 bits of x
 	target = ( argument.words.hi < signMask );         	// flag positive sign
-	
-	if ( xhi > 0x41e00000ul ) 
+
+	if ( xhi > 0x41e00000ul )
 /*******************************************************************************
 *     Is x is out of long range or NaN?                                        *
 *******************************************************************************/
@@ -327,19 +327,19 @@ long int roundtol ( double x )
 		else
 			return ( LONG_MIN );
 		}
-	
-	if ( target ) 
+
+	if ( target )
 /*******************************************************************************
 *     Is sign of x is "+"?                                                     *
 *******************************************************************************/
 		{
-		if ( x < 2147483647.5 ) 
+		if ( x < 2147483647.5 )
 /*******************************************************************************
 *     x is in the range of a long.                                             *
 *******************************************************************************/
 			{
 			y = ( x + doubleToLong ) - doubleToLong; 	// round at binary point
-			if ( y != x )	
+			if ( y != x )
 				{		                    	// inexact case
 				asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// save environment
 				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
@@ -348,7 +348,7 @@ long int roundtol ( double x )
 				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 				return ( ( long ) argument.words.lo );
 				}
-			
+
 			argument.dbl = y + doubleToLong;		// force result into argument.words.lo
 			return ( ( long ) argument.words.lo ); 	// return long result
 			}
@@ -363,13 +363,13 @@ long int roundtol ( double x )
 /*******************************************************************************
 *     x < 0.0 and may or may not be out of the range of a long.                *
 *******************************************************************************/
-	if ( x > -2147483648.5 ) 
+	if ( x > -2147483648.5 )
 /*******************************************************************************
 *     x is in the range of a long.                                             *
 *******************************************************************************/
 		{
 		y = ( x + doubleToLong ) - doubleToLong;	  	// round at binary point
-		if ( y != x ) 
+		if ( y != x )
 			{			                    	// inexact case
 			asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// save environment
 			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
@@ -378,7 +378,7 @@ long int roundtol ( double x )
 			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
 			return ( ( long ) argument.words.lo );
 			}
-		
+
 		argument.dbl = y + doubleToLong;
 		return ( ( long ) argument.words.lo );	  	//  return long result
 		}
@@ -398,29 +398,29 @@ long int roundtol ( double x )
 *     inexact if an ordered return value is not equal to the operand.          *
 *                                                                              *
 *******************************************************************************/
-   
+
 double trunc ( double x )
-      {	
+      {
 	DblInHex argument,OldEnvironment;
 	register double y;
 	register unsigned long int xhi;
 	register long int target;
-	
+
 	argument.dbl = x;
 	xhi = argument.words.hi & 0x7fffffffUL;	      	// xhi <- high half of |x|
 	target = ( argument.words.hi < signMask );	      	// flag positive sign
-	
-	if ( xhi < 0x43300000ul ) 
+
+	if ( xhi < 0x43300000ul )
 /*******************************************************************************
 *     Is |x| < 2.0^53?                                                         *
 *******************************************************************************/
 		{
-		if ( xhi < 0x3ff00000ul ) 
+		if ( xhi < 0x3ff00000ul )
 /*******************************************************************************
 *     Is |x| < 1.0?                                                            *
 *******************************************************************************/
 			{
-			if ( ( xhi | argument.words.lo ) != 0ul ) 
+			if ( ( xhi | argument.words.lo ) != 0ul )
 				{                             	// raise deserved INEXACT
 				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
 				OldEnvironment.words.lo |= 0x02000000ul;
@@ -434,7 +434,7 @@ double trunc ( double x )
 /*******************************************************************************
 *     Is 1.0 < |x| < 2.0^52?                                                   *
 *******************************************************************************/
-		if ( target ) 
+		if ( target )
 			{
 			y = ( x + twoTo52 ) - twoTo52;      	// round at binary point
 			if ( y > x )
@@ -442,8 +442,8 @@ double trunc ( double x )
 			else
 				return ( y );
 			}
-		
-		else 
+
+		else
 			{
 			y = ( x - twoTo52 ) + twoTo52;      	// round at binary point.
 			if ( y < x )
@@ -470,7 +470,7 @@ double trunc ( double x )
 *******************************************************************************/
 
 /*******************************************************************************
-*     modf is the double implementation.                                       *                             
+*     modf is the double implementation.                                       *
 *******************************************************************************/
 
 double modf ( double x, double *iptr )
@@ -478,19 +478,19 @@ double modf ( double x, double *iptr )
       register double OldEnvironment, xtrunc;
       register unsigned long int xHead, signBit;
       DblInHex argument;
-      
+
       argument.dbl = x;
       xHead = argument.words.hi & 0x7ffffffful;            // |x| high bit pattern
       signBit = ( argument.words.hi & 0x80000000ul );      // isolate sign bit
 	  if (xHead == 0x7ff81fe0)
 			signBit = signBit | 0;
-      
-      if ( xHead < 0x43300000ul ) 
+
+      if ( xHead < 0x43300000ul )
 /*******************************************************************************
 *     Is |x| < 2.0^53?                                                         *
 *******************************************************************************/
             {
-            if ( xHead < 0x3ff00000ul )      
+            if ( xHead < 0x3ff00000ul )
 /*******************************************************************************
 *     Is |x| < 1.0?                                                            *
 *******************************************************************************/
@@ -515,18 +515,18 @@ double modf ( double x, double *iptr )
             *iptr = xtrunc;                               // store integral part
             if ( x != xtrunc )                            // nonzero fraction
                   return ( x - xtrunc );
-            else 
+            else
                   {                                       // zero with x's sign
                   argument.words.hi = signBit;
                   argument.words.lo = 0ul;
                   return ( argument.dbl );
                   }
             }
-      
+
       *iptr = x;                                          // x is integral or NaN
       if ( x != x )                                       // NaN is returned
             return x;
-      else 
+      else
             {                                             // zero with x's sign
             argument.words.hi = signBit;
             argument.words.lo = 0ul;

libm/powerpc/s_rint.c

@@ -1,18 +1,18 @@
 /*******************************************************************************
 **      File:   rndint.c
-**      
+**
 **      Contains: C source code for implementations of floating-point
 **                functions which round to integral value or format, as
 **                defined in header <fp.h>.  In particular, this file
 **                contains implementations of functions rint, nearbyint,
 **                rinttol, round, roundtol, trunc, modf and modfl.  This file
 **                targets PowerPC or Power platforms.
-**                        
+**
 **      Written by: A. Sazegari, Apple AltiVec Group
 **	    Created originally by Jon Okada, Apple Numerics Group
-**      
+**
 **      Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
-**      
+**
 **      Change History (most recent first):
 **
 **      13 Jul 01  ram  replaced --setflm calls with inline assembly
@@ -21,7 +21,7 @@
 **				1. removed double_t, put in double for now.
 **				2. removed iclass from nearbyint.
 **				3. removed wrong comments intrunc.
-**				4. 
+**				4.
 **      13 May 97  ali  made performance improvements in rint, rinttol, roundtol
 **                      and trunc by folding some of the taligent ideas into this
 **                      implementation.  nearbyint is faster than the one in taligent,
@@ -33,7 +33,7 @@
 **      18 Feb 93  ali  Changed the return value of fenv functions
 **                      feclearexcept and feraiseexcept to their new
 **                      NCEG X3J11.1/93-001 definitions.
-**      16 Dec 92  JPO  Removed __itrunc implementation to a 
+**      16 Dec 92  JPO  Removed __itrunc implementation to a
 **                      separate file.
 **      15 Dec 92  JPO  Added __itrunc implementation and modified
 **                      rinttol to include conversion from double
@@ -41,7 +41,7 @@
 **                      call __itrunc.
 **      10 Dec 92  JPO  Added modf (double) implementation.
 **      04 Dec 92  JPO  First created.
-**                        
+**
 *******************************************************************************/
 
 #include <limits.h>
@@ -73,7 +73,7 @@ static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
 *                                                                              *
 *     The function rint rounds its double argument to integral value           *
 *     according to the current rounding direction and returns the result in    *
-*     double format.  This function signals inexact if an ordered return       * 
+*     double format.  This function signals inexact if an ordered return       *
 *     value is not equal to the operand.                                       *
 *                                                                              *
 ********************************************************************************
@@ -89,16 +89,16 @@ static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
 *double rint ( double x )
 *      {
 *      double y;
-*      
+*
 *      y = twoTo52.fval;
-*      
-*      if ( fabs ( x ) >= y )                          // huge case is exact 
+*
+*      if ( fabs ( x ) >= y )                          // huge case is exact
 *            return x;
-*      if ( x < 0 ) y = -y;                            // negative case 
-*      y = ( x + y ) - y;                              // force rounding 
-*      if ( y == 0.0 )                                 // zero results mirror sign of x 
+*      if ( x < 0 ) y = -y;                            // negative case
+*      y = ( x + y ) - y;                              // force rounding
+*      if ( y == 0.0 )                                 // zero results mirror sign of x
 *            y = copysign ( y, x );
-*      return ( y );      
+*      return ( y );
 *      }
 ********************************************************************************
 *     Now a bit twidling version that is about %30 faster.                     *
@@ -110,17 +110,17 @@ double rint ( double x )
       register double y;
       unsigned long int xHead;
       register long int target;
-      
+
       argument.dbl = x;
       xHead = argument.words.hi & 0x7fffffffUL;          // xHead <- high half of |x|
       target = ( argument.words.hi < signMask );         // flags positive sign
-      
-      if ( xHead < 0x43300000ul ) 
+
+      if ( xHead < 0x43300000ul )
 /*******************************************************************************
 *     Is |x| < 2.0^52?                                                         *
 *******************************************************************************/
             {
-            if ( xHead < 0x3ff00000ul ) 
+            if ( xHead < 0x3ff00000ul )
 /*******************************************************************************
 *     Is |x| < 1.0?                                                            *
 *******************************************************************************/
@@ -129,7 +129,7 @@ double rint ( double x )
                         y = ( x + twoTo52 ) - twoTo52;  // round at binary point
                   else
                         y = ( x - twoTo52 ) + twoTo52;  // round at binary point
-                  if ( y == 0.0 ) 
+                  if ( y == 0.0 )
                         {                               // fix sign of zero result
                         if ( target )
                               return ( 0.0 );
@@ -138,7 +138,7 @@ double rint ( double x )
                         }
                   return y;
                   }
-            
+
 /*******************************************************************************
 *     Is 1.0 < |x| < 2.0^52?                                                   *
 *******************************************************************************/
@@ -148,7 +148,7 @@ double rint ( double x )
             else
                   return ( ( x - twoTo52 ) + twoTo52 );
             }
-      
+
 /*******************************************************************************
 *     |x| >= 2.0^52 or x is a NaN.                                             *
 *******************************************************************************/

libm/powerpc/w_scalb.c

@@ -1,26 +1,26 @@
 /***********************************************************************
 **      File:    scalb.c
-**      
+**
 **      Contains: C source code for implementations of floating-point
 **                scalb functions defined in header <fp.h>.  In
 **                particular, this file contains implementations of
 **                functions scalb and scalbl for double and long double
 **                formats on PowerPC platforms.
-**                        
+**
 **      Written by: Jon Okada, SANEitation Engineer, ext. 4-4838
-**      
+**
 **      Copyright: © 1992 by Apple Computer, Inc., all rights reserved
-**      
+**
 **      Change History ( most recent first ):
 **
 **      28 May 97  ali   made an speed improvement for large n,
 **                       removed scalbl.
 **      12 Dec 92  JPO   First created.
-**                        
+**
 ***********************************************************************/
 
-typedef union           
-      { 
+typedef union
+      {
       struct {
 #if defined(__BIG_ENDIAN__)
         unsigned long int hi;
@@ -41,35 +41,35 @@ static const double twoToM1022 = 2.225073858507201383e-308;  // 0x1p-1022
       double  scalb( double  x, long int n ) returns its argument x scaled
       by the factor 2^m.  NaNs, signed zeros, and infinities are propagated
       by this function regardless of the value of n.
-      
+
       Exceptions:  OVERFLOW/INEXACT or UNDERFLOW inexact may occur;
                          INVALID for signaling NaN inputs ( quiet NaN returned ).
-      
+
       Calls:  none.
 ***********************************************************************/
 
 double scalb ( double x, int n  )
       {
       DblInHex xInHex;
-      
+
       xInHex.words.lo = 0UL;                     // init. low half of xInHex
-      
-      if ( n > 1023 ) 
+
+      if ( n > 1023 )
             {                                   // large positive scaling
             if ( n > 2097 )                     // huge scaling
             	return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023;
-            while ( n > 1023 ) 
+            while ( n > 1023 )
                   {                             // scale reduction loop
                   x *= twoTo1023;               // scale x by 2^1023
                   n -= 1023;                    // reduce n by 1023
                   }
             }
-      
-      else if ( n < -1022 ) 
+
+      else if ( n < -1022 )
             {                                   // large negative scaling
             if ( n < -2098 )                    // huge negative scaling
             	return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022;
-            while ( n < -1022 ) 
+            while ( n < -1022 )
                   {                             // scale reduction loop
                   x *= twoToM1022;              // scale x by 2^( -1022 )
                   n += 1022;                    // incr n by 1022

libm/s_asinh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -16,26 +16,26 @@ static char rcsid[] = "$NetBSD: s_asinh.c,v 1.9 1995/05/12 04:57:37 jtc Exp $";
 
 /* asinh(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)*(log(x)+ln2)) for large |x|, else
  *		 := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
- *		 := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))  
+ *		 := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
  */
 
 #include "math.h"
 #include "math_private.h"
 
 #ifdef __STDC__
-static const double 
+static const double
 #else
-static double 
+static double
 #endif
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
-huge=  1.00000000000000000000e+300; 
+huge=  1.00000000000000000000e+300;
 
 #ifdef __STDC__
 	double asinh(double x)
@@ -43,7 +43,7 @@ huge=  1.00000000000000000000e+300;
 	double asinh(x)
 	double x;
 #endif
-{	
+{
 	double t,w;
 	int32_t hx,ix;
 	GET_HIGH_WORD(hx,x);
@@ -51,7 +51,7 @@ huge=  1.00000000000000000000e+300;
 	if(ix>=0x7ff00000) return x+x;	/* x is inf or NaN */
 	if(ix< 0x3e300000) {	/* |x|<2**-28 */
 	    if(huge+x>one) return x;	/* return x inexact except 0 */
-	} 
+	}
 	if(ix>0x41b00000) {	/* |x| > 2**28 */
 	    w = __ieee754_log(fabs(x))+ln2;
 	} else if (ix>0x40000000) {	/* 2**28 > |x| > 2.0 */

libm/s_atan.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -28,9 +28,9 @@ static char rcsid[] = "$NetBSD: s_atan.c,v 1.8 1995/05/10 20:46:45 jtc Exp $";
  *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
- * compiler will convert from decimal to binary accurately enough 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
 
@@ -78,9 +78,9 @@ static double aT[] = {
 };
 
 #ifdef __STDC__
-	static const double 
+	static const double
 #else
-	static double 
+	static double
 #endif
 one   = 1.0,
 huge   = 1.0e300;
@@ -114,9 +114,9 @@ huge   = 1.0e300;
 	x = fabs(x);
 	if (ix < 0x3ff30000) {		/* |x| < 1.1875 */
 	    if (ix < 0x3fe60000) {	/* 7/16 <=|x|<11/16 */
-		id = 0; x = (2.0*x-one)/(2.0+x); 
+		id = 0; x = (2.0*x-one)/(2.0+x);
 	    } else {			/* 11/16<=|x|< 19/16 */
-		id = 1; x  = (x-one)/(x+one); 
+		id = 1; x  = (x-one)/(x+one);
 	    }
 	} else {
 	    if (ix < 0x40038000) {	/* |x| < 2.4375 */

libm/s_cbrt.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -40,9 +40,9 @@ F =  1.60714285714285720630e+00, /* 45/28     = 0x3FF9B6DB, 0x6DB6DB6E */
 G =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
 
 #ifdef __STDC__
-	double cbrt(double x) 
+	double cbrt(double x)
 #else
-	double cbrt(x) 
+	double cbrt(x)
 	double x;
 #endif
 {
@@ -56,7 +56,7 @@ G =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
 	hx  ^=sign;
 	if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
 	GET_LOW_WORD(low,x);
-	if((hx|low)==0) 
+	if((hx|low)==0)
 	    return(x);		/* cbrt(0) is itself */
 
 	SET_HIGH_WORD(x,hx);	/* x <- |x| */
@@ -72,9 +72,9 @@ G =  3.57142857142857150787e-01; /* 5/14      = 0x3FD6DB6D, 0xB6DB6DB7 */
     /* 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) */
 	GET_HIGH_WORD(high,t);
 	INSERT_WORDS(t,high+0x00000001,0);
 

libm/s_ceil.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -46,7 +46,7 @@ static double huge = 1.0e300;
 	if(j0<20) {
 	    if(j0<0) { 	/* raise inexact if x != 0 */
 		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
-		    if(i0<0) {i0=0x80000000;i1=0;} 
+		    if(i0<0) {i0=0x80000000;i1=0;}
 		    else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
 		}
 	    } else {
@@ -65,7 +65,7 @@ static double huge = 1.0e300;
 	    if((i1&i)==0) return x;	/* x is integral */
 	    if(huge+x>0.0) { 		/* raise inexact flag */
 		if(i0>0) {
-		    if(j0==20) i0+=1; 
+		    if(j0==20) i0+=1;
 		    else {
 			j = i1 + (1<<(52-j0));
 			if(j<i1) i0+=1;	/* got a carry */

libm/s_ceilf.c

@@ -8,7 +8,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -41,7 +41,7 @@ static float huge = 1.0e30;
 	if(j0<23) {
 	    if(j0<0) { 	/* raise inexact if x != 0 */
 		if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
-		    if(i0<0) {i0=0x80000000;} 
+		    if(i0<0) {i0=0x80000000;}
 		    else if(i0!=0) { i0=0x3f800000;}
 		}
 	    } else {

libm/s_copysign.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/s_cos.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -23,8 +23,8 @@ static char rcsid[] = "$NetBSD: s_cos.c,v 1.7 1995/05/10 20:47:02 jtc Exp $";
  *	__ieee754_rem_pio2	... argument reduction routine
  *
  * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on 
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
  *	in [-pi/4 , +pi/4], and let n = k mod 4.
  *	We have
  *
@@ -42,7 +42,7 @@ static char rcsid[] = "$NetBSD: s_cos.c,v 1.7 1995/05/10 20:47:02 jtc Exp $";
  *      trig(NaN)    is that NaN;
  *
  * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded 
+ *	TRIG(x) returns trig(x) nearly rounded
  */
 
 #include "math.h"

libm/s_erf.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -19,11 +19,11 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
  *			     x
  *		      2      |\
  *     erf(x)  =  ---------  | exp(-t*t)dt
- *	 	   sqrt(pi) \| 
+ *	 	   sqrt(pi) \|
  *			     0
  *
  *     erfc(x) =  1-erf(x)
- *  Note that 
+ *  Note that
  *		erf(-x) = -erf(x)
  *		erfc(-x) = 2 - erfc(x)
  *
@@ -36,7 +36,7 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
  *	   Q is an odd poly of degree 10.
  *						 -57.90
  *			| R - (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 + ....)
@@ -59,14 +59,14 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
  *	   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
  *
- *      3. For x in [1.25,1/0.35(~2.857143)], 
+ *      3. For x in [1.25,1/0.35(~2.857143)],
  *         	erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
  *         	erf(x)  = 1 - erfc(x)
- *	   where 
+ *	   where
  *		R1(z) = degree 7 poly in z, (z=1/x^2)
  *		S1(z) = degree 8 poly in z
  *
@@ -84,7 +84,7 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
  *	   To compute exp(-x*x-0.5625+R/S), let s be a single
  *	   precision number and s := x; then
  *		-x*x = -s*s + (s-x)*(s+x)
- *	        exp(-x*x-0.5626+R/S) = 
+ *	        exp(-x*x-0.5626+R/S) =
  *			exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
  *      Note2:
  *	   Here 4 and 5 make use of the asymptotic series
@@ -104,7 +104,7 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $";
  *
  *      7. Special case:
  *         	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
  */
 
@@ -139,7 +139,7 @@ qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
 qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
 qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
 /*
- * Coefficients for approximation to  erf  in [0.84375,1.25] 
+ * Coefficients for approximation to  erf  in [0.84375,1.25]
  */
 pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
 pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
@@ -192,9 +192,9 @@ sb6  =  4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
 sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 
 #ifdef __STDC__
-	double erf(double x) 
+	double erf(double x)
 #else
-	double erf(x) 
+	double erf(x)
 	double x;
 #endif
 {
@@ -209,7 +209,7 @@ sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 
 	if(ix < 0x3feb0000) {		/* |x|<0.84375 */
 	    if(ix < 0x3e300000) { 	/* |x|<2**-28 */
-	        if (ix < 0x00800000) 
+	        if (ix < 0x00800000)
 		    return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
 		return x + efx*x;
 	    }
@@ -241,16 +241,16 @@ sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
 				sb5+s*(sb6+s*sb7))))));
 	}
-	z  = x;  
+	z  = x;
 	SET_LOW_WORD(z,0);
 	r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
 	if(hx>=0) return one-r/x; else return  r/x-one;
 }
 
 #ifdef __STDC__
-	double erfc(double x) 
+	double erfc(double x)
 #else
-	double erfc(x) 
+	double erfc(x)
 	double x;
 #endif
 {
@@ -283,7 +283,7 @@ sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 	    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(hx>=0) {
-	        z  = one-erx; return z - P/Q; 
+	        z  = one-erx; return z - P/Q;
 	    } else {
 		z = erx+P/Q; return one+z;
 	    }

libm/s_expm1.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -21,9 +21,9 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $";
  *   1. Argument reduction:
  *	Given x, find r and integer k such that
  *
- *               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658  
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658
  *
- *      Here a correction term c will be computed to compensate 
+ *      Here a correction term c will be computed to compensate
  *	the error in r when rounded to a floating-point number.
  *
  *   2. Approximating expm1(r) by a special rational function on
@@ -36,9 +36,9 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $";
  *	    R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
  *		     = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
  *		     = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
- *      We use a special Reme algorithm on [0,0.347] to generate 
- * 	a polynomial of degree 5 in r*r to approximate R1. The 
- *	maximum error of this polynomial approximation is bounded 
+ *      We use a special Reme algorithm on [0,0.347] to generate
+ * 	a polynomial of degree 5 in r*r to approximate R1. The
+ *	maximum error of this polynomial approximation is bounded
  *	by 2**-61. In other words,
  *	    R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
  *	where 	Q1  =  -1.6666666666666567384E-2,
@@ -49,28 +49,28 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $";
  *  	(where z=r*r, and the values of Q1 to Q5 are listed below)
  *	with error bounded by
  *	    |                  5           |     -61
- *	    | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2 
+ *	    | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
  *	    |                              |
- *	
- *	expm1(r) = exp(r)-1 is then computed by the following 
- * 	specific way which minimize the accumulation rounding error: 
+ *
+ *	expm1(r) = exp(r)-1 is then computed by the following
+ * 	specific way which minimize the accumulation rounding error:
  *			       2     3
  *			      r     r    [ 3 - (R1 + R1*r/2)  ]
  *	      expm1(r) = r + --- + --- * [--------------------]
  *		              2     2    [ 6 - r*(3 - R1*r/2) ]
- *	
+ *
  *	To compensate the error in the argument reduction, we use
- *		expm1(r+c) = expm1(r) + c + expm1(r)*c 
- *			   ~ expm1(r) + c + r*c 
+ *		expm1(r+c) = expm1(r) + c + expm1(r)*c
+ *			   ~ expm1(r) + c + r*c
  *	Thus c+r*c will be added in as the correction terms for
- *	expm1(r+c). Now rearrange the term to avoid optimization 
+ *	expm1(r+c). Now rearrange the term to avoid optimization
  * 	screw up:
  *		        (      2                                    2 )
  *		        ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
  *	 expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
  *	                ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
  *                      (                                             )
- *    	
+ *
  *		   = r - E
  *   3. Scale back to obtain expm1(x):
  *	From step 1, we have
@@ -87,7 +87,7 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $";
  *	       	       else	     return  1.0+2.0*(r-E);
  *	  (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
  *	  (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
- *	  (vii) return 2^k(1-((E+2^-k)-r)) 
+ *	  (vii) return 2^k(1-((E+2^-k)-r))
  *
  * Special cases:
  *	expm1(INF) is INF, expm1(NaN) is NaN;
@@ -99,12 +99,12 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $";
  *	1 ulp (unit in the last place).
  *
  * Misc. info.
- *	For IEEE double 
+ *	For IEEE double
  *	    if x >  7.09782712893383973096e+02 then expm1(x) overflow
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
  * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  */
@@ -153,7 +153,7 @@ Q5  =  -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
                 if(hx>=0x7ff00000) {
 		    u_int32_t low;
 		    GET_LOW_WORD(low,x);
-		    if(((hx&0xfffff)|low)!=0) 
+		    if(((hx&0xfffff)|low)!=0)
 		         return x+x; 	 /* NaN */
 		    else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
 	        }
@@ -166,7 +166,7 @@ Q5  =  -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
 	}
 
     /* argument reduction */
-	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */ 
+	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */
 	    if(hx < 0x3FF0A2B2) {	/* and |x| < 1.5 ln2 */
 		if(xsb==0)
 		    {hi = x - ln2_hi; lo =  ln2_lo;  k =  1;}
@@ -180,10 +180,10 @@ Q5  =  -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
 	    }
 	    x  = hi - lo;
 	    c  = (hi-x)-lo;
-	} 
+	}
 	else if(hx < 0x3c900000) {  	/* when |x|<2**-54, return x */
 	    t = huge+x;	/* return x with inexact flags when x!=0 */
-	    return x - (t-(huge+x));	
+	    return x - (t-(huge+x));
 	}
 	else k = 0;
 
@@ -198,7 +198,7 @@ Q5  =  -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
 	    e  = (x*(e-c)-c);
 	    e -= hxs;
 	    if(k== -1) return 0.5*(x-e)-0.5;
-	    if(k==1) { 
+	    if(k==1) {
 	       	if(x < -0.25) return -2.0*(e-(x+0.5));
 	       	else 	      return  one+2.0*(x-e);
 	    }

libm/s_fabs.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/s_finite.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/s_floor.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -46,7 +46,7 @@ static double huge = 1.0e300;
 	if(j0<20) {
 	    if(j0<0) { 	/* raise inexact if x != 0 */
 		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
-		    if(i0>=0) {i0=i1=0;} 
+		    if(i0>=0) {i0=i1=0;}
 		    else if(((i0&0x7fffffff)|i1)!=0)
 			{ i0=0xbff00000;i1=0;}
 		}
@@ -66,7 +66,7 @@ static double huge = 1.0e300;
 	    if((i1&i)==0) return x;	/* x is integral */
 	    if(huge+x>0.0) { 		/* raise inexact flag */
 		if(i0<0) {
-		    if(j0==20) i0+=1; 
+		    if(j0==20) i0+=1;
 		    else {
 			j = i1+(1<<(52-j0));
 			if(j<i1) i0 +=1 ; 	/* got a carry */

libm/s_floorf.c

@@ -8,7 +8,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -49,7 +49,7 @@ static float huge = 1.0e30;
 	if(j0<23) {
 	    if(j0<0) { 	/* raise inexact if x != 0 */
 		if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
-		    if(i0>=0) {i0=0;} 
+		    if(i0>=0) {i0=0;}
 		    else if((i0&0x7fffffff)!=0)
 			{ i0=0xbf800000;}
 		}

libm/s_frexp.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,13 +15,13 @@ static char rcsid[] = "$NetBSD: s_frexp.c,v 1.9 1995/05/10 20:47:24 jtc Exp $";
 #endif
 
 /*
- * for non-zero x 
+ * for non-zero x
  *	x = frexp(arg,&exp);
  * return a double fp quantity x such that 0.5 <= |x| <1.0
  * and the corresponding binary exponent "exp". That is
  *	arg = x*2^exp.
- * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg 
- * with *exp=0. 
+ * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
+ * with *exp=0.
  */
 
 #include "math.h"

libm/s_ilogb.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -36,7 +36,7 @@ static char rcsid[] = "$NetBSD: s_ilogb.c,v 1.9 1995/05/10 20:47:28 jtc Exp $";
 	hx &= 0x7fffffff;
 	if(hx<0x00100000) {
 	    GET_LOW_WORD(lx,x);
-	    if((hx|lx)==0) 
+	    if((hx|lx)==0)
 		return 0x80000001;	/* ilogb(0) = 0x80000001 */
 	    else			/* subnormal x */
 		if(hx==0) {

libm/s_ldexp.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/s_lib_version.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/s_log1p.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -16,9 +16,9 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $";
 
 /* double log1p(double x)
  *
- * Method :                  
- *   1. Argument Reduction: find k and f such that 
- *			1+x = 2^k * (1+f), 
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *			1+x = 2^k * (1+f),
  *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
  *
  *      Note. If k=0, then f=x is exact. However, if k!=0, then f
@@ -32,8 +32,8 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $";
  *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
  *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
  *	     	 = 2s + s*R
- *      We use a special Reme algorithm on [0,0.1716] to generate 
- * 	a polynomial of degree 14 to approximate R The maximum error 
+ *      We use a special Reme algorithm on [0,0.1716] to generate
+ * 	a polynomial of degree 14 to approximate R The maximum error
  *	of this polynomial approximation is bounded by 2**-58.45. In
  *	other words,
  *		        2      4      6      8      10      12      14
@@ -41,21 +41,21 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $";
  *  	(the values of Lp1 to Lp7 are listed in the program)
  *	and
  *	    |      2          14          |     -58.45
- *	    | Lp1*s +...+Lp7*s    -  R(z) | <= 2 
+ *	    | Lp1*s +...+Lp7*s    -  R(z) | <= 2
  *	    |                             |
  *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
  *	In order to guarantee error in log below 1ulp, we compute log
  *	by
  *		log1p(f) = f - (hfsq - s*(hfsq+R)).
- *	
- *	3. Finally, log1p(x) = k*ln2 + log1p(f).  
+ *
+ *	3. Finally, log1p(x) = k*ln2 + log1p(f).
  *		 	     = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
- *	   Here ln2 is split into two floating point number: 
+ *	   Here ln2 is split into two floating point number:
  *			ln2_hi + ln2_lo,
  *	   where n*ln2_hi is always exact for |n| < 2000.
  *
  * Special cases:
- *	log1p(x) is NaN with signal if x < -1 (including -INF) ; 
+ *	log1p(x) is NaN with signal if x < -1 (including -INF) ;
  *	log1p(+INF) is +INF; log1p(-1) is -INF with signal;
  *	log1p(NaN) is that NaN with no signal.
  *
@@ -64,14 +64,14 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $";
  *	1 ulp (unit in the last place).
  *
  * Constants:
- * The hexadecimal values are the intended ones for the following 
- * constants. The decimal values may be used, provided that the 
- * compiler will convert from decimal to binary accurately enough 
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
  * to produce the hexadecimal values shown.
  *
  * Note: Assuming log() return accurate answer, the following
  * 	 algorithm can be used to compute log1p(x) to within a few ULP:
- *	
+ *
  *		u = 1+x;
  *		if(u==1.0) return x ; else
  *			   return log(u)*(x/(u-1.0));
@@ -132,11 +132,11 @@ static double zero = 0.0;
 	    }
 	    if(hx>0||hx<=((int32_t)0xbfd2bec3)) {
 		k=0;f=x;hu=1;}	/* -0.2929<x<0.41422 */
-	} 
+	}
 	if (hx >= 0x7ff00000) return x+x;
 	if(k!=0) {
 	    if(hx<0x43400000) {
-		u  = 1.0+x; 
+		u  = 1.0+x;
 		GET_HIGH_WORD(hu,u);
 	        k  = (hu>>20)-1023;
 	        c  = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
@@ -151,7 +151,7 @@ static double zero = 0.0;
 	    if(hu<0x6a09e) {
 	        SET_HIGH_WORD(u,hu|0x3ff00000);	/* normalize u */
 	    } else {
-	        k += 1; 
+	        k += 1;
 		SET_HIGH_WORD(u,hu|0x3fe00000);	/* normalize u/2 */
 	        hu = (0x00100000-hu)>>2;
 	    }
@@ -159,14 +159,14 @@ static double zero = 0.0;
 	}
 	hfsq=0.5*f*f;
 	if(hu==0) {	/* |f| < 2**-20 */
-	    if(f==zero) {if(k==0) return zero;  
+	    if(f==zero) {if(k==0) return zero;
 			else {c += k*ln2_lo; return k*ln2_hi+c;}
 	    }
 	    R = hfsq*(1.0-0.66666666666666666*f);
 	    if(k==0) return f-R; else
 	    	     return k*ln2_hi-((R-(k*ln2_lo+c))-f);
 	}
- 	s = f/(2.0+f); 
+ 	s = f/(2.0+f);
 	z = s*s;
 	R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
 	if(k==0) return f-(hfsq-s*(hfsq+R)); else

libm/s_logb.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -36,7 +36,7 @@ static char rcsid[] = "$NetBSD: s_logb.c,v 1.8 1995/05/10 20:47:50 jtc Exp $";
 	if((ix|lx)==0) return -1.0/fabs(x);
 	if(ix>=0x7ff00000) return x*x;
 	if((ix>>=20)==0) 			/* IEEE 754 logb */
-		return -1022.0; 
+		return -1022.0;
 	else
-		return (double) (ix-1023); 
+		return (double) (ix-1023);
 }

libm/s_matherr.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/s_modf.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -15,7 +15,7 @@ static char rcsid[] = "$NetBSD: s_modf.c,v 1.8 1995/05/10 20:47:55 jtc Exp $";
 #endif
 
 /*
- * modf(double x, double *iptr) 
+ * modf(double x, double *iptr)
  * return fraction part of x, and return x's integral part in *iptr.
  * Method:
  *	Bit twiddling.

libm/s_nextafter.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -39,15 +39,15 @@ static char rcsid[] = "$NetBSD: s_nextafter.c,v 1.8 1995/05/10 20:47:58 jtc Exp
 	ix = hx&0x7fffffff;		/* |x| */
 	iy = hy&0x7fffffff;		/* |y| */
 
-	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */ 
-	   ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0))     /* y is nan */ 
-	   return x+y;				
+	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */
+	   ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0))     /* y is nan */
+	   return x+y;
 	if(x==y) return x;		/* x=y, return x */
 	if((ix|lx)==0) {			/* x == 0 */
 	    INSERT_WORDS(x,hy&0x80000000,1);	/* return +-minsubnormal */
 	    y = x*x;
 	    if(y==x) return y; else return x;	/* raise underflow flag */
-	} 
+	}
 	if(hx>=0) {				/* x > 0 */
 	    if(hx>hy||((hx==hy)&&(lx>ly))) {	/* x > y, x -= ulp */
 		if(lx==0) hx -= 1;

libm/s_rint.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -30,7 +30,7 @@ static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
 #ifdef __STDC__
 static const double
 #else
-static double 
+static double
 #endif
 TWO52[2]={
   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
@@ -51,7 +51,7 @@ TWO52[2]={
 	sx = (i0>>31)&1;
 	j0 = ((i0>>20)&0x7ff)-0x3ff;
 	if(j0<20) {
-	    if(j0<0) { 	
+	    if(j0<0) {
 		if(((i0&0x7fffffff)|i1)==0) return x;
 		i1 |= (i0&0x0fffff);
 		i0 &= 0xfffe0000;

libm/s_scalbn.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,10 +14,10 @@
 static char rcsid[] = "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
 #endif
 
-/* 
+/*
  * scalbn (double x, int n)
- * scalbn(x,n) returns x* 2**n  computed by  exponent  
- * manipulation rather than by actually performing an 
+ * scalbn(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
  * exponentiation or a multiplication.
  */
 
@@ -46,13 +46,13 @@ tiny   = 1.0e-300;
         k = (hx&0x7ff00000)>>20;		/* extract exponent */
         if (k==0) {				/* 0 or subnormal x */
             if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
-	    x *= two54; 
+	    x *= two54;
 	    GET_HIGH_WORD(hx,x);
-	    k = ((hx&0x7ff00000)>>20) - 54; 
+	    k = ((hx&0x7ff00000)>>20) - 54;
             if (n< -50000) return tiny*x; 	/*underflow*/
 	    }
         if (k==0x7ff) return x+x;		/* NaN or Inf */
-        k = k+n; 
+        k = k+n;
         if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
         if (k > 0) 				/* normal result */
 	    {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}

libm/s_significand.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/s_sin.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -23,8 +23,8 @@ static char rcsid[] = "$NetBSD: s_sin.c,v 1.7 1995/05/10 20:48:15 jtc Exp $";
  *	__ieee754_rem_pio2	... argument reduction routine
  *
  * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on 
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
  *	in [-pi/4 , +pi/4], and let n = k mod 4.
  *	We have
  *
@@ -42,7 +42,7 @@ static char rcsid[] = "$NetBSD: s_sin.c,v 1.7 1995/05/10 20:48:15 jtc Exp $";
  *      trig(NaN)    is that NaN;
  *
  * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded 
+ *	TRIG(x) returns trig(x) nearly rounded
  */
 
 #include "math.h"

libm/s_tan.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -22,8 +22,8 @@ static char rcsid[] = "$NetBSD: s_tan.c,v 1.7 1995/05/10 20:48:18 jtc Exp $";
  *	__ieee754_rem_pio2	... argument reduction routine
  *
  * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on 
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
+ *      Let S,C and T denote the sin, cos and tan respectively on
+ *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
  *	in [-pi/4 , +pi/4], and let n = k mod 4.
  *	We have
  *
@@ -41,7 +41,7 @@ static char rcsid[] = "$NetBSD: s_tan.c,v 1.7 1995/05/10 20:48:18 jtc Exp $";
  *      trig(NaN)    is that NaN;
  *
  * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded 
+ *	TRIG(x) returns trig(x) nearly rounded
  */
 
 #include "math.h"

libm/s_tanh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -62,7 +62,7 @@ static double one=1.0, two=2.0, tiny = 1.0e-300;
 	ix = jx&0x7fffffff;
 
     /* x is INF or NaN */
-	if(ix>=0x7ff00000) { 
+	if(ix>=0x7ff00000) {
 	    if (jx>=0) return one/x+one;    /* tanh(+-inf)=+-1 */
 	    else       return one/x-one;    /* tanh(NaN) = NaN */
 	}

libm/w_acos.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/w_acosh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_acosh.c,v 1.6 1995/05/10 20:48:31 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper acosh(x)
  */
 

libm/w_asin.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_asin.c,v 1.6 1995/05/10 20:48:35 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper asin(x)
  */
 

libm/w_atan2.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_atan2.c,v 1.6 1995/05/10 20:48:39 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper atan2(y,x)
  */
 #include "math.h"

libm/w_atanh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_atanh.c,v 1.6 1995/05/10 20:48:43 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper atanh(x)
  */
 
@@ -39,7 +39,7 @@ static char rcsid[] = "$NetBSD: w_atanh.c,v 1.6 1995/05/10 20:48:43 jtc Exp $";
 	if(y>=1.0) {
 	    if(y>1.0)
 	        return __kernel_standard(x,x,30); /* atanh(|x|>1) */
-	    else 
+	    else
 	        return __kernel_standard(x,x,31); /* atanh(|x|==1) */
 	} else
 	    return z;

libm/w_cabs.c

@@ -1,6 +1,6 @@
 /*
  * cabs() wrapper for hypot().
- * 
+ *
  * Written by J.T. Conklin, <jtc@wimsey.com>
  * Placed into the Public Domain, 1994.
  */

libm/w_cosh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_cosh.c,v 1.6 1995/05/10 20:48:47 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper cosh(x)
  */
 
@@ -34,7 +34,7 @@ static char rcsid[] = "$NetBSD: w_cosh.c,v 1.6 1995/05/10 20:48:47 jtc Exp $";
 	double z;
 	z = __ieee754_cosh(x);
 	if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
-	if(fabs(x)>7.10475860073943863426e+02) {	
+	if(fabs(x)>7.10475860073943863426e+02) {
 	        return __kernel_standard(x,x,5); /* cosh overflow */
 	} else
 	    return z;

libm/w_drem.c

@@ -1,6 +1,6 @@
 /*
  * drem() wrapper for remainder().
- * 
+ *
  * Written by J.T. Conklin, <jtc@wimsey.com>
  * Placed into the Public Domain, 1994.
  */

libm/w_exp.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_exp.c,v 1.6 1995/05/10 20:48:51 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper exp(x)
  */
 
@@ -47,7 +47,7 @@ u_threshold= -7.45133219101941108420e+02;  /* 0xc0874910, 0xD52D3051 */
 	        return __kernel_standard(x,x,6); /* exp overflow */
 	    else if(x<u_threshold)
 	        return __kernel_standard(x,x,7); /* exp underflow */
-	} 
+	}
 	return z;
 #endif
 }

libm/w_fmod.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_fmod.c,v 1.6 1995/05/10 20:48:55 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper fmod(x,y)
  */
 

libm/w_gamma.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -46,4 +46,4 @@ extern int signgam;
         } else
             return y;
 #endif
-}             
+}

libm/w_gamma_r.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_gamma_r.c,v 1.7 1995/11/20 22:06:45 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper double gamma_r(double x, int *signgamp)
  */
 
@@ -43,4 +43,4 @@ static char rcsid[] = "$NetBSD: w_gamma_r.c,v 1.7 1995/11/20 22:06:45 jtc Exp $"
         } else
             return y;
 #endif
-}             
+}

libm/w_hypot.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/w_j0.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */

libm/w_j1.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,8 +14,8 @@
 static char rcsid[] = "$NetBSD: w_j1.c,v 1.6 1995/05/10 20:49:15 jtc Exp $";
 #endif
 
-/* 
- * wrapper of j1,y1 
+/*
+ * wrapper of j1,y1
  */
 
 #include "math.h"

libm/w_jn.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -18,7 +18,7 @@ static char rcsid[] = "$NetBSD: w_jn.c,v 1.6 1995/05/10 20:49:19 jtc Exp $";
  * wrapper 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.
@@ -37,7 +37,7 @@ static char rcsid[] = "$NetBSD: w_jn.c,v 1.6 1995/05/10 20:49:19 jtc Exp $";
  *	yn(n,x) is similar in all respects, except
  *	that forward recursion is used for all
  *	values of n>1.
- *	
+ *
  */
 
 #include "math.h"

libm/w_lgamma.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -46,4 +46,4 @@ extern int signgam;
         } else
             return y;
 #endif
-}             
+}

libm/w_lgamma_r.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_lgamma_r.c,v 1.6 1995/05/10 20:49:27 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper double lgamma_r(double x, int *signgamp)
  */
 
@@ -43,4 +43,4 @@ static char rcsid[] = "$NetBSD: w_lgamma_r.c,v 1.6 1995/05/10 20:49:27 jtc Exp $
         } else
             return y;
 #endif
-}             
+}

libm/w_log.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -37,7 +37,7 @@ static char rcsid[] = "$NetBSD: w_log.c,v 1.6 1995/05/10 20:49:33 jtc Exp $";
 	if(_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0) return z;
 	if(x==0.0)
 	    return __kernel_standard(x,x,16); /* log(0) */
-	else 
+	else
 	    return __kernel_standard(x,x,17); /* log(x<0) */
 #endif
 }

libm/w_log10.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_log10.c,v 1.6 1995/05/10 20:49:35 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper log10(X)
  */
 
@@ -38,7 +38,7 @@ static char rcsid[] = "$NetBSD: w_log10.c,v 1.6 1995/05/10 20:49:35 jtc Exp $";
 	if(x<=0.0) {
 	    if(x==0.0)
 	        return __kernel_standard(x,x,18); /* log10(0) */
-	    else 
+	    else
 	        return __kernel_standard(x,x,19); /* log10(x<0) */
 	} else
 	    return z;

libm/w_pow.c

@@ -7,12 +7,12 @@
  *
  * 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.
  * ====================================================
  */
 
-/* 
+/*
  * wrapper pow(x,y) return x**y
  */
 
@@ -34,12 +34,12 @@
 	z=__ieee754_pow(x,y);
 	if(_LIB_VERSION == _IEEE_|| isnan(y)) return z;
 	if(isnan(x)) {
-	    if(y==0.0) 
+	    if(y==0.0)
 	        return __kernel_standard(x,y,42); /* pow(NaN,0.0) */
-	    else 
+	    else
 		return z;
 	}
-	if(x==0.0){ 
+	if(x==0.0){
 	    if(y==0.0)
 	        return __kernel_standard(x,y,20); /* pow(0.0,0.0) */
 	    if(finite(y)&&y<0.0)
@@ -50,10 +50,10 @@
 	    if(finite(x)&&finite(y)) {
 	        if(isnan(z))
 	            return __kernel_standard(x,y,24); /* pow neg**non-int */
-	        else 
+	        else
 	            return __kernel_standard(x,y,21); /* pow overflow */
 	    }
-	} 
+	}
 	if(z==0.0&&finite(x)&&finite(y))
 	    return __kernel_standard(x,y,22); /* pow underflow */
 	return z;

libm/w_remainder.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_remainder.c,v 1.6 1995/05/10 20:49:44 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper remainder(x,p)
  */
 
@@ -34,7 +34,7 @@ static char rcsid[] = "$NetBSD: w_remainder.c,v 1.6 1995/05/10 20:49:44 jtc Exp
 	double z;
 	z = __ieee754_remainder(x,y);
 	if(_LIB_VERSION == _IEEE_ || isnan(y)) return z;
-	if(y==0.0) 
+	if(y==0.0)
 	    return __kernel_standard(x,y,28); /* remainder(x,0) */
 	else
 	    return z;

libm/w_scalb.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -16,7 +16,7 @@ static char rcsid[] = "$NetBSD: w_scalb.c,v 1.6 1995/05/10 20:49:48 jtc Exp $";
 
 /*
  * wrapper scalb(double x, double fn) is provide for
- * passing various standard test suite. One 
+ * passing various standard test suite. One
  * should use scalbn() instead.
  */
 
@@ -51,10 +51,10 @@ static char rcsid[] = "$NetBSD: w_scalb.c,v 1.6 1995/05/10 20:49:48 jtc Exp $";
 	}
 	if(z==0.0&&z!=x) {
 	    return __kernel_standard(x,(double)fn,33); /* scalb underflow */
-	} 
+	}
 #ifndef _SCALB_INT
 	if(!finite(fn)) errno = ERANGE;
 #endif
 	return z;
-#endif 
+#endif
 }

libm/w_sinh.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_sinh.c,v 1.6 1995/05/10 20:49:51 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper sinh(x)
  */
 
@@ -31,7 +31,7 @@ static char rcsid[] = "$NetBSD: w_sinh.c,v 1.6 1995/05/10 20:49:51 jtc Exp $";
 #ifdef _IEEE_LIBM
 	return __ieee754_sinh(x);
 #else
-	double z; 
+	double z;
 	z = __ieee754_sinh(x);
 	if(_LIB_VERSION == _IEEE_) return z;
 	if(!finite(z)&&finite(x)) {

libm/w_sqrt.c

@@ -5,7 +5,7 @@
  *
  * Developed at SunPro, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
  */
@@ -14,7 +14,7 @@
 static char rcsid[] = "$NetBSD: w_sqrt.c,v 1.6 1995/05/10 20:49:55 jtc Exp $";
 #endif
 
-/* 
+/*
  * wrapper sqrt(x)
  */
 

libm/w_sqrtf.c

@@ -18,7 +18,7 @@
  * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  */
 
-/* 
+/*
  * wrapper for sqrt(x)
  */