diff options
Diffstat (limited to 'lib/libm/common_source')
| -rw-r--r-- | lib/libm/common_source/acosh.c | 10 | ||||
| -rw-r--r-- | lib/libm/common_source/asincos.c | 38 | ||||
| -rw-r--r-- | lib/libm/common_source/asinh.c | 12 | ||||
| -rw-r--r-- | lib/libm/common_source/atan.c | 18 | ||||
| -rw-r--r-- | lib/libm/common_source/atanh.c | 4 | ||||
| -rw-r--r-- | lib/libm/common_source/cosh.c | 26 | ||||
| -rw-r--r-- | lib/libm/common_source/erf.c | 16 | ||||
| -rw-r--r-- | lib/libm/common_source/exp.c | 20 | ||||
| -rw-r--r-- | lib/libm/common_source/exp__E.c | 10 | ||||
| -rw-r--r-- | lib/libm/common_source/expm1.c | 32 | ||||
| -rw-r--r-- | lib/libm/common_source/j0.c | 28 | ||||
| -rw-r--r-- | lib/libm/common_source/j1.c | 32 | ||||
| -rw-r--r-- | lib/libm/common_source/jn.c | 44 | ||||
| -rw-r--r-- | lib/libm/common_source/log.c | 2 | ||||
| -rw-r--r-- | lib/libm/common_source/log10.c | 10 | ||||
| -rw-r--r-- | lib/libm/common_source/log1p.c | 30 | ||||
| -rw-r--r-- | lib/libm/common_source/log__L.c | 12 | ||||
| -rw-r--r-- | lib/libm/common_source/pow.c | 20 | ||||
| -rw-r--r-- | lib/libm/common_source/sinh.c | 8 | ||||
| -rw-r--r-- | lib/libm/common_source/tanh.c | 4 | 
20 files changed, 188 insertions, 188 deletions
| diff --git a/lib/libm/common_source/acosh.c b/lib/libm/common_source/acosh.c index bc16cc7b46a8..149e5deec1e0 100644 --- a/lib/libm/common_source/acosh.c +++ b/lib/libm/common_source/acosh.c @@ -48,10 +48,10 @@ static char sccsid[] = "@(#)acosh.c	8.1 (Berkeley) 6/4/93";   *	log1p(x) 		...return log(1+x)   *   * Method : - *	Based on  + *	Based on   *		acosh(x) = log [ x + sqrt(x*x-1) ]   *	we have - *		acosh(x) := log1p(x)+ln2,	if (x > 1.0E20); else		 + *		acosh(x) := log1p(x)+ln2,	if (x > 1.0E20); else   *		acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) .   *	These formulae avoid the over/underflow complication.   * @@ -60,7 +60,7 @@ static char sccsid[] = "@(#)acosh.c	8.1 (Berkeley) 6/4/93";   *	acosh(NaN) is NaN without signal.   *   * Accuracy: - *	acosh(x) returns the exact inverse hyperbolic cosine of x nearly  + *	acosh(x) returns the exact inverse hyperbolic cosine of x nearly   *	rounded. In a test run with 512,000 random arguments on a VAX, the   *	maximum observed error was 3.30 ulps (units of the last place) at   *	x=1.0070493753568216 . @@ -87,7 +87,7 @@ ic(ln2lo, 1.9082149292705877000E-10,-33, 1.A39EF35793C76)  double acosh(x)  double x; -{	 +{  	double t,big=1.E20; /* big+1==big */  #if !defined(vax)&&!defined(tahoe) @@ -95,7 +95,7 @@ double x;  #endif	/* !defined(vax)&&!defined(tahoe) */      /* return log1p(x) + log(2) if x is large */ -	if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);}  +	if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);}  	t=sqrt(x-1.0);  	return(log1p(t*(t+sqrt(x+1.0)))); diff --git a/lib/libm/common_source/asincos.c b/lib/libm/common_source/asincos.c index c746b1652ac2..12d0e14e1ae3 100644 --- a/lib/libm/common_source/asincos.c +++ b/lib/libm/common_source/asincos.c @@ -45,12 +45,12 @@ static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";   *	sqrt(x)   *   * Required kernel function: - *	atan2(y,x)  + *	atan2(y,x)   * - * Method :                   - *	asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is  + * Method : + *	asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is   *		  computed as follows - *			1-x*x                     if x <  0.5,  + *			1-x*x                     if x <  0.5,   *			2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5.   *   * Special cases: @@ -59,22 +59,22 @@ static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";   *   * Accuracy:   * 1)  If atan2() uses machine PI, then - *  + *   *	asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded;   *	and PI is the exact pi rounded to machine precision (see atan2 for   *      details):   *   *	in decimal: - *		pi = 3.141592653589793 23846264338327 .....  + *		pi = 3.141592653589793 23846264338327 .....   *    53 bits   PI = 3.141592653589793 115997963 ..... , - *    56 bits   PI = 3.141592653589793 227020265 ..... ,   + *    56 bits   PI = 3.141592653589793 227020265 ..... ,   *   *	in hexadecimal:   *		pi = 3.243F6A8885A308D313198A2E....   *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps   *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps - *	 - *	In a test run with more than 200,000 random arguments on a VAX, the  + * + *	In a test run with more than 200,000 random arguments on a VAX, the   *	maximum observed error in ulps (units in the last place) was   *	2.06 ulps.      (comparing against (PI/pi)*(exact asin(x)));   * @@ -82,7 +82,7 @@ static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";   *   *	asin(x) returns the exact asin(x) with error below about 2 ulps.   * - *	In a test run with more than 1,024,000 random arguments on a VAX, the  + *	In a test run with more than 1,024,000 random arguments on a VAX, the   *	maximum observed error in ulps (units in the last place) was   *      1.99 ulps.   */ @@ -97,7 +97,7 @@ double x;  	s=copysign(x,one);  	if(s <= 0.5)  	    return(atan2(x,sqrt(one-x*x))); -	else  +	else  	    { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); }  } @@ -112,9 +112,9 @@ double x;   *	sqrt(x)   *   * Required kernel function: - *	atan2(y,x)  + *	atan2(y,x)   * - * Method :                   + * Method :   *			      ________   *                           / 1 - x   *	acos(x) = 2*atan2(  / -------- , 1 ) . @@ -126,22 +126,22 @@ double x;   *   * Accuracy:   * 1)  If atan2() uses machine PI, then - *  + *   *	acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded;   *	and PI is the exact pi rounded to machine precision (see atan2 for   *      details):   *   *	in decimal: - *		pi = 3.141592653589793 23846264338327 .....  + *		pi = 3.141592653589793 23846264338327 .....   *    53 bits   PI = 3.141592653589793 115997963 ..... , - *    56 bits   PI = 3.141592653589793 227020265 ..... ,   + *    56 bits   PI = 3.141592653589793 227020265 ..... ,   *   *	in hexadecimal:   *		pi = 3.243F6A8885A308D313198A2E....   *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps   *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps - *	 - *	In a test run with more than 200,000 random arguments on a VAX, the  + * + *	In a test run with more than 200,000 random arguments on a VAX, the   *	maximum observed error in ulps (units in the last place) was   *	2.07 ulps.      (comparing against (PI/pi)*(exact acos(x)));   * @@ -149,7 +149,7 @@ double x;   *   *	acos(x) returns the exact acos(x) with error below about 2 ulps.   * - *	In a test run with more than 1,024,000 random arguments on a VAX, the  + *	In a test run with more than 1,024,000 random arguments on a VAX, the   *	maximum observed error in ulps (units in the last place) was   *	2.15 ulps.   */ diff --git a/lib/libm/common_source/asinh.c b/lib/libm/common_source/asinh.c index 5db8d2ddf7f8..1804145e3a93 100644 --- a/lib/libm/common_source/asinh.c +++ b/lib/libm/common_source/asinh.c @@ -49,16 +49,16 @@ static char sccsid[] = "@(#)asinh.c	8.1 (Berkeley) 6/4/93";   *	log1p(x) 		...return log(1+x)   *   * Method : - *	Based on  + *	Based on   *		asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]   *	we have   *	asinh(x) := x  if  1+x*x=1,   *		 := sign(x)*(log1p(x)+ln2))	 if sqrt(1+x*x)=x, else - *		 := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) )   + *		 := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) )   *   * Accuracy:   *	asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded. - *	In a test run with 52,000 random arguments on a VAX, the maximum  + *	In a test run with 52,000 random arguments on a VAX, the maximum   *	observed error was 1.58 ulps (units in the last place).   *   * Constants: @@ -82,16 +82,16 @@ ic(ln2lo, 1.9082149292705877000E-10, -33, 1.A39EF35793C76)  double asinh(x)  double x; -{	 +{  	double t,s;  	const static double	small=1.0E-10,	/* fl(1+small*small) == 1 */  				big  =1.0E20,	/* fl(1+big) == big */ -				one  =1.0   ;	 +				one  =1.0   ;  #if !defined(vax)&&!defined(tahoe)  	if(x!=x) return(x);	/* x is NaN */  #endif	/* !defined(vax)&&!defined(tahoe) */ -	if((t=copysign(x,one))>small)  +	if((t=copysign(x,one))>small)  	    if(t<big) {  	     	s=one/t; return(copysign(log1p(t+t/(s+sqrt(one+s*s))),x)); }  	    else	/* if |x| > big */ diff --git a/lib/libm/common_source/atan.c b/lib/libm/common_source/atan.c index 272c7f1cdc48..f29c7d4d3483 100644 --- a/lib/libm/common_source/atan.c +++ b/lib/libm/common_source/atan.c @@ -41,32 +41,32 @@ static char sccsid[] = "@(#)atan.c	8.1 (Berkeley) 6/4/93";   * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85.   *   * Required kernel function: - *	atan2(y,x)  + *	atan2(y,x)   * - * Method:                   - *	atan(x) = atan2(x,1.0).  + * Method: + *	atan(x) = atan2(x,1.0).   *   * Special case:   *	if x is NaN, return x itself.   *   * Accuracy:   * 1)  If atan2() uses machine PI, then - *  + *   *	atan(x) returns (PI/pi) * (the exact arc tangent of x) nearly rounded;   *	and PI is the exact pi rounded to machine precision (see atan2 for   *      details):   *   *	in decimal: - *		pi = 3.141592653589793 23846264338327 .....  + *		pi = 3.141592653589793 23846264338327 .....   *    53 bits   PI = 3.141592653589793 115997963 ..... , - *    56 bits   PI = 3.141592653589793 227020265 ..... ,   + *    56 bits   PI = 3.141592653589793 227020265 ..... ,   *   *	in hexadecimal:   *		pi = 3.243F6A8885A308D313198A2E....   *    53 bits   PI = 3.243F6A8885A30  =  2 * 1.921FB54442D18	error=.276ulps   *    56 bits   PI = 3.243F6A8885A308 =  4 * .C90FDAA22168C2    error=.206ulps - *	 - *	In a test run with more than 200,000 random arguments on a VAX, the  + * + *	In a test run with more than 200,000 random arguments on a VAX, the   *	maximum observed error in ulps (units in the last place) was   *	0.86 ulps.      (comparing against (PI/pi)*(exact atan(x))).   * @@ -74,7 +74,7 @@ static char sccsid[] = "@(#)atan.c	8.1 (Berkeley) 6/4/93";   *   *	atan(x) returns the exact atan(x) with error below about 2 ulps.   * - *	In a test run with more than 1,024,000 random arguments on a VAX, the  + *	In a test run with more than 1,024,000 random arguments on a VAX, the   *	maximum observed error in ulps (units in the last place) was   *	0.85 ulps.   */ diff --git a/lib/libm/common_source/atanh.c b/lib/libm/common_source/atanh.c index 89cb61cca25e..e5cdadde9052 100644 --- a/lib/libm/common_source/atanh.c +++ b/lib/libm/common_source/atanh.c @@ -38,14 +38,14 @@ static char sccsid[] = "@(#)atanh.c	8.1 (Berkeley) 6/4/93";  /* ATANH(X)   * RETURN THE HYPERBOLIC ARC TANGENT OF X   * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85;  + * CODED IN C BY K.C. NG, 1/8/85;   * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85.   *   * Required kernel function:   *	log1p(x) 	...return log(1+x)   *   * Method : - *	Return  + *	Return   *                          1              2x                          x   *		atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)   *                          2             1 - x                      1 - x diff --git a/lib/libm/common_source/cosh.c b/lib/libm/common_source/cosh.c index e2b30731b8c1..e8d35195eb05 100644 --- a/lib/libm/common_source/cosh.c +++ b/lib/libm/common_source/cosh.c @@ -38,7 +38,7 @@ static char sccsid[] = "@(#)cosh.c	8.1 (Berkeley) 6/4/93";  /* COSH(X)   * RETURN THE HYPERBOLIC COSINE OF X   * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85;  + * CODED IN C BY K.C. NG, 1/8/85;   * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85.   *   * Required system supported functions : @@ -46,20 +46,20 @@ static char sccsid[] = "@(#)cosh.c	8.1 (Berkeley) 6/4/93";   *	scalb(x,N)   *   * Required kernel function: - *	exp(x)  + *	exp(x)   *	exp__E(x,c)	...return exp(x+c)-1-x for |x|<0.3465   *   * Method : - *	1. Replace x by |x|.  - *	2.  - *		                                        [ exp(x) - 1 ]^2  + *	1. Replace x by |x|. + *	2. + *		                                        [ exp(x) - 1 ]^2   *	    0        <= x <= 0.3465  :  cosh(x) := 1 + -------------------   *			       			           2*exp(x)   *   *		                                   exp(x) +  1/exp(x)   *	    0.3465   <= x <= 22      :  cosh(x) := -------------------   *			       			           2 - *	    22       <= x <= lnovfl  :  cosh(x) := exp(x)/2  + *	    22       <= x <= lnovfl  :  cosh(x) := exp(x)/2   *	    lnovfl   <= x <= lnovfl+log(2)   *				     :  cosh(x) := exp(x)/2 (avoid overflow)   *	    log(2)+lnovfl <  x <  INF:  overflow to INF @@ -106,7 +106,7 @@ static max = 1023                     ;  double cosh(x)  double x; -{	 +{  	static const double half=1.0/2.0,  		one=1.0, small=1.0E-18; /* fl(1+small)==1 */  	double t; @@ -115,19 +115,19 @@ double x;  	if(x!=x) return(x);	/* x is NaN */  #endif	/* !defined(vax)&&!defined(tahoe) */  	if((x=copysign(x,one)) <= 22) -	    if(x<0.3465)  +	    if(x<0.3465)  		if(x<small) return(one+x);  		else {t=x+__exp__E(x,0.0);x=t+t; return(one+t*t/(2.0+x)); }  	    else /* for x lies in [0.3465,22] */  	        { t=exp(x); return((t+one/t)*half); } -	if( lnovfl <= x && x <= (lnovfl+0.7))  -        /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1))  -         * and return 2^max*exp(x) to avoid unnecessary overflow  +	if( lnovfl <= x && x <= (lnovfl+0.7)) +        /* for x lies in [lnovfl, lnovfl+ln2], decrease x by ln(2^(max+1)) +         * and return 2^max*exp(x) to avoid unnecessary overflow           */ -	    return(scalb(exp((x-mln2hi)-mln2lo), max));  +	    return(scalb(exp((x-mln2hi)-mln2lo), max)); -	else  +	else  	    return(exp(x)*half);	/* for large x,  cosh(x)=exp(x)/2 */  } diff --git a/lib/libm/common_source/erf.c b/lib/libm/common_source/erf.c index 308f1a952d73..ba0c006eea5f 100644 --- a/lib/libm/common_source/erf.c +++ b/lib/libm/common_source/erf.c @@ -68,13 +68,13 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";   *          erfc(x) = 1 - erf(x)           if x<=0.25   *                  = 0.5 + ((0.5-x)-x*P)  if x in [0.25,0.84375]   *	   where - *			2		 2	  4		  20   + *			2		 2	  4		  20   *              P =  P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x  )   * 	   is an approximation to (erf(x)-x)/x with precision   *   *						 -56.45   *			| P - (erf(x)-x)/x | <= 2 - *	 + *   *   *	   Remark. The formula is derived by noting   *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) @@ -96,7 +96,7 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";   *	   That is, we use rational approximation to approximate   *			erf(1+s) - (c = (single)0.84506291151)   *	   Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] - *	   where  + *	   where   *		P1(s) = degree 6 poly in s   *		Q1(s) = degree 6 poly in s   * @@ -105,7 +105,7 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";   *		erfc(x)	= (1/x)exp(-x*x-(.5*log(pi) -.5z + R(z)/S(z))   *   *	Where z = 1/(x*x), R is degree 9, and S is degree 3; - *	 + *   *      5. For x in [4,28]   *         	erf(x)  = 1.0 - tiny   *		erfc(x)	= (1/x)exp(-x*x-(.5*log(pi)+eps + zP(z)) @@ -139,7 +139,7 @@ static char sccsid[] = "@(#)erf.c	8.1 (Berkeley) 6/4/93";   *   *      7. Special cases:   *         	erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1, - *         	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,  + *         	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,   *	   	erfc/erf(NaN) is NaN   */ @@ -181,7 +181,7 @@ p8 =   1.640186161764254363152286358441771740838e-0006,  p9 =  -1.571599331700515057841960987689515895479e-0007,  p10=   1.073087585213621540635426191486561494058e-0008;  /* - * Coefficients for approximation to erf in [0.84375,1.25]  + * Coefficients for approximation to erf in [0.84375,1.25]   */  static double  pa0 =  -2.362118560752659485957248365514511540287e-0003, @@ -319,7 +319,7 @@ double erf(x)  		return (z-one);  } -double erfc(x)  +double erfc(x)  	double x;  {  	double R,S,P,Q,s,ax,y,z,r,fabs(),__exp__D(); @@ -352,7 +352,7 @@ double erfc(x)  	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));  	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));  	    if (x>=0) { -	        z  = one-c; return z - P/Q;  +	        z  = one-c; return z - P/Q;  	    } else {  		z = c+P/Q; return one+z;  	    } diff --git a/lib/libm/common_source/exp.c b/lib/libm/common_source/exp.c index 9b4f045f82e4..fa6ea75735c5 100644 --- a/lib/libm/common_source/exp.c +++ b/lib/libm/common_source/exp.c @@ -38,21 +38,21 @@ static char sccsid[] = "@(#)exp.c	8.1 (Berkeley) 6/4/93";  /* EXP(X)   * RETURN THE EXPONENTIAL OF X   * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) - * CODED IN C BY K.C. NG, 1/19/85;  + * CODED IN C BY K.C. NG, 1/19/85;   * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.   *   * Required system supported functions: - *	scalb(x,n)	 - *	copysign(x,y)	 + *	scalb(x,n) + *	copysign(x,y)   *	finite(x)   *   * Method: - *	1. Argument Reduction: given the input x, find r and integer k such  + *	1. Argument Reduction: given the input x, find r and integer k such   *	   that - *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .   + *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .   *	   r will be represented as r := z+c for better accuracy.   * - *	2. Compute exp(r) by  + *	2. Compute exp(r) by   *   *		exp(r) = 1 + r + r*R1/(2-R1),   *	   where @@ -143,7 +143,7 @@ double x;  		}  		/* end of x > lntiny */ -		else  +		else  		     /* exp(-big#) underflows to zero */  		     if(finite(x))  return(scalb(1.0,-5000)); @@ -152,7 +152,7 @@ double x;  	}  	/* end of x < lnhuge */ -	else  +	else  	/* exp(INF) is INF, exp(+big#) overflows to INF */  	    return( finite(x) ?  scalb(1.0,5000)  : x);  } @@ -188,7 +188,7 @@ double x, c;  		}  		/* end of x > lntiny */ -		else  +		else  		     /* exp(-big#) underflows to zero */  		     if(finite(x))  return(scalb(1.0,-5000)); @@ -197,7 +197,7 @@ double x, c;  	}  	/* end of x < lnhuge */ -	else  +	else  	/* exp(INF) is INF, exp(+big#) overflows to INF */  	    return( finite(x) ?  scalb(1.0,5000)  : x);  } diff --git a/lib/libm/common_source/exp__E.c b/lib/libm/common_source/exp__E.c index ab972477a04b..16ec4cbb2179 100644 --- a/lib/libm/common_source/exp__E.c +++ b/lib/libm/common_source/exp__E.c @@ -41,7 +41,7 @@ static char sccsid[] = "@(#)exp__E.c	8.1 (Berkeley) 6/4/93";   * exp__E RETURNS   *   *			 /  exp(x+c) - 1 - x ,  1E-19 < |x| < .3465736 - *       exp__E(x,c) = 	| 		      + *       exp__E(x,c) = 	|   *			 \  0 ,  |x| < 1E-19.   *   * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) @@ -50,12 +50,12 @@ static char sccsid[] = "@(#)exp__E.c	8.1 (Berkeley) 6/4/93";   * REVISED BY K.C. NG on 3/16/85, 4/16/85.   *   * Required system supported function: - *	copysign(x,y)	 + *	copysign(x,y)   *   * Method:   *	1. Rational approximation. Let r=x+c.   *	   Based on - *                                   2 * sinh(r/2)      + *                                   2 * sinh(r/2)   *                exp(r) - 1 =   ----------------------   ,   *                               cosh(r/2) - sinh(r/2)   *	   exp__E(r) is computed using @@ -76,7 +76,7 @@ static char sccsid[] = "@(#)exp__E.c	8.1 (Berkeley) 6/4/93";   * Approximation error:   *   *   |	exp(x) - 1			   |        2**(-57),  (IEEE double) - *   | ------------  -  (exp__E(x,0)+x)/x  |  <=  + *   | ------------  -  (exp__E(x,0)+x)/x  |  <=   *   |	     x			           |	    2**(-69).  (VAX D)   *   * Constants: @@ -120,7 +120,7 @@ double x,c;  #else	/* defined(vax)||defined(tahoe) */             q = z*( q1 +z*  q2 );  #endif	/* defined(vax)||defined(tahoe) */ -           xp= x*p     ;  +           xp= x*p     ;  	   xh= x*half  ;             w = xh-(q-xp)  ;  	   p = p+p; diff --git a/lib/libm/common_source/expm1.c b/lib/libm/common_source/expm1.c index 760d2bed677d..8a84f14c4da2 100644 --- a/lib/libm/common_source/expm1.c +++ b/lib/libm/common_source/expm1.c @@ -38,36 +38,36 @@ static char sccsid[] = "@(#)expm1.c	8.1 (Berkeley) 6/4/93";  /* EXPM1(X)   * RETURN THE EXPONENTIAL OF X MINUS ONE   * DOUBLE PRECISION (IEEE 53 BITS, VAX D FORMAT 56 BITS) - * CODED IN C BY K.C. NG, 1/19/85;  + * CODED IN C BY K.C. NG, 1/19/85;   * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/21/85, 4/16/85.   *   * Required system supported functions: - *	scalb(x,n)	 - *	copysign(x,y)	 + *	scalb(x,n) + *	copysign(x,y)   *	finite(x)   *   * Kernel function:   *	exp__E(x,c)   *   * Method: - *	1. Argument Reduction: given the input x, find r and integer k such  + *	1. Argument Reduction: given the input x, find r and integer k such   *	   that - *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .   + *	                   x = k*ln2 + r,  |r| <= 0.5*ln2 .   *	   r will be represented as r := z+c for better accuracy.   * - *	2. Compute EXPM1(r)=exp(r)-1 by  + *	2. Compute EXPM1(r)=exp(r)-1 by   *   *			EXPM1(r=z+c) := z + exp__E(z,c)   *   *	3. EXPM1(x) =  2^k * ( EXPM1(r) + 1-2^-k ).   * - * 	Remarks:  + * 	Remarks:   *	   1. When k=1 and z < -0.25, we use the following formula for   *	      better accuracy:   *			EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) )   *	   2. To avoid rounding error in 1-2^-k where k is large, we use   *			EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 } - *	      when k>56.  + *	      when k>56.   *   * Special cases:   *	EXPM1(INF) is INF, EXPM1(NaN) is NaN; @@ -108,7 +108,7 @@ ic(invln2, 1.4426950408889633870E0,     0, 1.71547652B82FE)  double expm1(x)  double x;  { -	const static double one=1.0, half=1.0/2.0;  +	const static double one=1.0, half=1.0/2.0;  	double  z,hi,lo,c;  	int k;  #if defined(vax)||defined(tahoe) @@ -126,13 +126,13 @@ double x;  		    /* argument reduction : x - k*ln2 */  			k= invln2 *x+copysign(0.5,x);	/* k=NINT(x/ln2) */ -			hi=x-k*ln2hi ;  +			hi=x-k*ln2hi ;  			z=hi-(lo=k*ln2lo);  			c=(hi-z)-lo;  			if(k==0) return(z+__exp__E(z,c));  			if(k==1) -			    if(z< -0.25)  +			    if(z< -0.25)  				{x=z+half;x +=__exp__E(z,c); return(x+x);}  			    else  				{z+=__exp__E(z,c); x=half+z; return(x+x);} @@ -143,17 +143,17 @@ double x;  			      { x=one-scalb(one,-k); z += __exp__E(z,c);}  			    else if(k<100)  			      { x = __exp__E(z,c)-scalb(one,-k); x+=z; z=one;} -			    else  +			    else  			      { x = __exp__E(z,c)+z; z=one;} -			    return (scalb(x+z,k));   +			    return (scalb(x+z,k));  			}  		}  		/* end of x > lnunfl */ -		else  +		else  		     /* expm1(-big#) rounded to -1 (inexact) */ -		     if(finite(x))   +		     if(finite(x))  			 { ln2hi+ln2lo; return(-one);}  		     /* expm1(-INF) is -1 */ @@ -161,7 +161,7 @@ double x;  	}  	/* end of x < lnhuge */ -	else  +	else  	/*  expm1(INF) is INF, expm1(+big#) overflows to INF */  	    return( finite(x) ?  scalb(one,5000) : x);  } diff --git a/lib/libm/common_source/j0.c b/lib/libm/common_source/j0.c index a9b28b38caea..684fb4d42c3e 100644 --- a/lib/libm/common_source/j0.c +++ b/lib/libm/common_source/j0.c @@ -46,18 +46,18 @@ static char sccsid[] = "@(#)j0.c	8.2 (Berkeley) 11/30/93";   *   * Developed at SunPro, a Sun Microsystems, Inc. business.   * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice  + * software is freely granted, provided that this notice   * is preserved.   * ====================================================   *   * ******************* WARNING ********************   * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision  + * Distributable Math Library) for IEEE double precision   * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE  - * Standard 754/854. This alpha version is distributed  - * for testing purpose. Those who use this software  - * should report any bugs to  + * in C that runs on machines that conform to IEEE + * Standard 754/854. This alpha version is distributed + * for testing purpose. Those who use this software + * should report any bugs to   *   *		fdlibm-comments@sunpro.eng.sun.com   * @@ -84,20 +84,20 @@ static char sccsid[] = "@(#)j0.c	8.2 (Berkeley) 11/30/93";   * 	   (To avoid cancellation, use   *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))   * 	    to compute the worse one.) - *	    + *   *	3 Special cases   *		j0(nan)= nan   *		j0(0) = 1   *		j0(inf) = 0 - *		 + *   * Method -- y0(x):   *	1. For x<2. - *	   Since  + *	   Since   *		y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)   *	   therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.   *	   We use the following function to approximate y0,   *		y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 - *	   where  + *	   where   *		U(z) = u0 + u1*z + ... + u6*z^6   *		V(z) = 1  + v1*z + ... + v4*z^4   *	   with absolute approximation error bounded by 2**-72. @@ -121,7 +121,7 @@ static char sccsid[] = "@(#)j0.c	8.2 (Berkeley) 11/30/93";  static double pzero __P((double)), qzero __P((double)); -static double  +static double  huge 	= 1e300,  zero    = 0.0,  one	= 1.0, @@ -138,7 +138,7 @@ s03 =   5.135465502073181376284426245689510134134e-0007,  s04 =   1.166140033337900097836930825478674320464e-0009;  double -j0(x)  +j0(x)  	double x;  {  	double z, s,c,ss,cc,r,u,v; @@ -201,7 +201,7 @@ v03 =   2.591508518404578033173189144579208685163e-0007,  v04 =   4.411103113326754838596529339004302243157e-0010;  double -y0(x)  +y0(x)  	double x;  {  	double z, s, c, ss, cc, u, v; @@ -345,7 +345,7 @@ static double pzero(x)  	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));  	return one+ r/s;  } -		 +  /* For x >= 8, the asymptotic expansions of qzero is   *	-1/8 s + 75/1024 s^3 - ..., where s = 1/x. diff --git a/lib/libm/common_source/j1.c b/lib/libm/common_source/j1.c index 71602aac1381..e8ca43ad83b8 100644 --- a/lib/libm/common_source/j1.c +++ b/lib/libm/common_source/j1.c @@ -46,18 +46,18 @@ static char sccsid[] = "@(#)j1.c	8.2 (Berkeley) 11/30/93";   *   * Developed at SunPro, a Sun Microsystems, Inc. business.   * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice  + * software is freely granted, provided that this notice   * is preserved.   * ====================================================   *   * ******************* WARNING ********************   * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision  + * Distributable Math Library) for IEEE double precision   * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE  - * Standard 754/854. This alpha version is distributed  - * for testing purpose. Those who use this software  - * should report any bugs to  + * in C that runs on machines that conform to IEEE + * Standard 754/854. This alpha version is distributed + * for testing purpose. Those who use this software + * should report any bugs to   *   *		fdlibm-comments@sunpro.eng.sun.com   * @@ -85,16 +85,16 @@ static char sccsid[] = "@(#)j1.c	8.2 (Berkeley) 11/30/93";   * 	   (To avoid cancellation, use   *		sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))   * 	    to compute the worse one.) - *	    + *   *	3 Special cases   *		j1(nan)= nan   *		j1(0) = 0   *		j1(inf) = 0 - *		 + *   * Method -- y1(x): - *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN  + *	1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN   *	2. For x<2. - *	   Since  + *	   Since   *		y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)   *	   therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.   *	   We use the following function to approximate y1, @@ -122,7 +122,7 @@ static char sccsid[] = "@(#)j1.c	8.2 (Berkeley) 11/30/93";  static double pone(), qone(); -static double  +static double  huge    = 1e300,  zero    = 0.0,  one	= 1.0, @@ -142,7 +142,7 @@ s05 =   1.235422744261379203512624973117299248281e-0011;  #define two_129	6.80564733841876926e+038	/* 2^129 */  #define two_m54	5.55111512312578270e-017	/* 2^-54 */ -double j1(x)  +double j1(x)  	double x;  {  	double z, s,c,ss,cc,r,u,v,y; @@ -205,7 +205,7 @@ static double v0[5] = {     1.665592462079920695971450872592458916421e-0011,  }; -double y1(x)  +double y1(x)  	double x;  {  	double z, s, c, ss, cc, u, v; @@ -254,10 +254,10 @@ double y1(x)                      z = invsqrtpi*(u*ss+v*cc)/sqrt(x);                  }                  return z; -        }  +        }          if (x <= two_m54) {    /* x < 2**-54 */              return (-tpi/x); -        }  +        }          z = x*x;          u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4])));          v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])))); @@ -351,7 +351,7 @@ static double pone(x)  	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));  	return (one + r/s);  } -		 +  /* For x >= 8, the asymptotic expansions of qone is   *	3/8 s - 105/1024 s^3 - ..., where s = 1/x. diff --git a/lib/libm/common_source/jn.c b/lib/libm/common_source/jn.c index 85a54012ecae..28d9687a51b1 100644 --- a/lib/libm/common_source/jn.c +++ b/lib/libm/common_source/jn.c @@ -46,18 +46,18 @@ static char sccsid[] = "@(#)jn.c	8.2 (Berkeley) 11/30/93";   *   * Developed at SunPro, a Sun Microsystems, Inc. business.   * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice  + * software is freely granted, provided that this notice   * is preserved.   * ====================================================   *   * ******************* WARNING ********************   * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision  + * Distributable Math Library) for IEEE double precision   * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE  - * Standard 754/854. This alpha version is distributed  - * for testing purpose. Those who use this software  - * should report any bugs to  + * in C that runs on machines that conform to IEEE + * Standard 754/854. This alpha version is distributed + * for testing purpose. Those who use this software + * should report any bugs to   *   *		fdlibm-comments@sunpro.eng.sun.com   * @@ -69,7 +69,7 @@ static char sccsid[] = "@(#)jn.c	8.2 (Berkeley) 11/30/93";   * jn(int n, double x), yn(int n, double x)   * floating point Bessel's function of the 1st and 2nd kind   * of order n - *           + *   * Special cases:   *	y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;   *	y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. @@ -88,7 +88,7 @@ static char sccsid[] = "@(#)jn.c	8.2 (Berkeley) 11/30/93";   *	yn(n,x) is similar in all respects, except   *	that forward recursion is used for all   *	values of n>1. - *	 + *   */  #include <math.h> @@ -120,7 +120,7 @@ double jn(n,x)       */      /* if J(n,NaN) is NaN */  	if (_IEEE && isnan(x)) return x+x; -	if (n<0){		 +	if (n<0){  		n = -n;  		x = -x;  	} @@ -134,10 +134,10 @@ double jn(n,x)  			/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */  	    if (_IEEE && x >= 8.148143905337944345e+090) {  					/* x >= 2**302 */ -    /* (x >> n**2)  +    /* (x >> n**2)       *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)       *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) -     *	    Let s=sin(x), c=cos(x),  +     *	    Let s=sin(x), c=cos(x),       *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then       *       *		   n	sin(xn)*sqt2	cos(xn)*sqt2 @@ -154,7 +154,7 @@ double jn(n,x)  		    case 3: temp =  cos(x)-sin(x); break;  		}  		b = invsqrtpi*temp/sqrt(x); -	    } else {	 +	    } else {  	        a = j0(x);  	        b = j1(x);  	        for(i=1;i<n;i++){ @@ -165,7 +165,7 @@ double jn(n,x)  	    }  	} else {  	    if (x < 1.86264514923095703125e-009) { /* x < 2**-29 */ -    /* x is tiny, return the first Taylor expansion of J(n,x)  +    /* x is tiny, return the first Taylor expansion of J(n,x)       * J(n,x) = 1/n!*(x/2)^n  - ...       */  		if (n > 33)	/* underflow */ @@ -180,14 +180,14 @@ double jn(n,x)  		}  	    } else {  		/* use backward recurrence */ -		/* 			x      x^2      x^2        +		/* 			x      x^2      x^2  		 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....  		 *			2n  - 2(n+1) - 2(n+2)  		 * -		 * 			1      1        1        +		 * 			1      1        1  		 *  (for large x)   =  ----  ------   ------   .....  		 *			2n   2(n+1)   2(n+2) -		 *			-- - ------ - ------ -  +		 *			-- - ------ - ------ -  		 *			 x     x         x  		 *  		 * Let w = 2n/x and h=2/x, then the above quotient @@ -203,9 +203,9 @@ double jn(n,x)  		 * To determine how many terms needed, let  		 * Q(0) = w, Q(1) = w(w+h) - 1,  		 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), -		 * When Q(k) > 1e4	good for single  -		 * When Q(k) > 1e9	good for double  -		 * When Q(k) > 1e17	good for quadruple  +		 * When Q(k) > 1e4	good for single +		 * When Q(k) > 1e9	good for double +		 * When Q(k) > 1e17	good for quadruple  		 */  	    /* determine k */  		double t,v; @@ -254,7 +254,7 @@ double jn(n,x)  	}  	return ((sgn == 1) ? -b : b);  } -double yn(n,x)  +double yn(n,x)  	int n; double x;  {  	int i, sign; @@ -275,10 +275,10 @@ double yn(n,x)  	if (n == 0) return(y0(x));  	if (n == 1) return(sign*y1(x));  	if(_IEEE && x >= 8.148143905337944345e+090) { /* x > 2**302 */ -    /* (x >> n**2)  +    /* (x >> n**2)       *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)       *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) -     *	    Let s=sin(x), c=cos(x),  +     *	    Let s=sin(x), c=cos(x),       *		xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then       *       *		   n	sin(xn)*sqt2	cos(xn)*sqt2 diff --git a/lib/libm/common_source/log.c b/lib/libm/common_source/log.c index ae186722f8df..908b8544f93c 100644 --- a/lib/libm/common_source/log.c +++ b/lib/libm/common_source/log.c @@ -391,7 +391,7 @@ log(x) double x;  			return (x+x);  		else  			return (infnan(ERANGE)); -	 +  	/* Argument reduction: 1 <= g < 2; x/2^m = g;	*/  	/* y = F*(1 + f/F) for |f| <= 2^-8		*/ diff --git a/lib/libm/common_source/log10.c b/lib/libm/common_source/log10.c index d2617ccebb9e..75205cd5ee15 100644 --- a/lib/libm/common_source/log10.c +++ b/lib/libm/common_source/log10.c @@ -38,9 +38,9 @@ static char sccsid[] = "@(#)log10.c	8.1 (Berkeley) 6/4/93";  /* LOG10(X)   * RETURN THE BASE 10 LOGARITHM OF x   * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/20/85;  + * CODED IN C BY K.C. NG, 1/20/85;   * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85. - *  + *   * Required kernel function:   *	log(x)   * @@ -52,12 +52,12 @@ static char sccsid[] = "@(#)log10.c	8.1 (Berkeley) 6/4/93";   *    Note:   *	  [log(10)]   rounded to 56 bits has error  .0895  ulps,   *	  [1/log(10)] rounded to 53 bits has error  .198   ulps; - *	  therefore, for better accuracy, in VAX D format, we divide  - *	  log(x) by log(10), but in IEEE Double format, we multiply  + *	  therefore, for better accuracy, in VAX D format, we divide + *	  log(x) by log(10), but in IEEE Double format, we multiply   *	  log(x) by [1/log(10)].   *   * Special cases: - *	log10(x) is NaN with signal if x < 0;  + *	log10(x) is NaN with signal if x < 0;   *	log10(+INF) is +INF with no signal; log10(0) is -INF with signal;   *	log10(NaN) is that NaN with no signal.   * diff --git a/lib/libm/common_source/log1p.c b/lib/libm/common_source/log1p.c index cbf9fcd895e2..020266721ed6 100644 --- a/lib/libm/common_source/log1p.c +++ b/lib/libm/common_source/log1p.c @@ -35,24 +35,24 @@  static char sccsid[] = "@(#)log1p.c	8.1 (Berkeley) 6/4/93";  #endif /* not lint */ -/* LOG1P(x)  +/* LOG1P(x)   * RETURN THE LOGARITHM OF 1+x   * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/19/85;  + * CODED IN C BY K.C. NG, 1/19/85;   * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85. - *  + *   * Required system supported functions: - *	scalb(x,n)  + *	scalb(x,n)   *	copysign(x,y) - *	logb(x)	 + *	logb(x)   *	finite(x)   *   * Required kernel function:   *	log__L(z)   *   * Method : - *	1. Argument Reduction: find k and f such that  - *			1+x  = 2^k * (1+f),  + *	1. Argument Reduction: find k and f such that + *			1+x  = 2^k * (1+f),   *	   where  sqrt(2)/2 < 1+f < sqrt(2) .   *   *	2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) @@ -65,11 +65,11 @@ static char sccsid[] = "@(#)log1p.c	8.1 (Berkeley) 6/4/93";   *   *	   See log__L() for the values of the coefficients.   * - *	3. Finally,  log(1+x) = k*ln2 + log(1+f).   + *	3. Finally,  log(1+x) = k*ln2 + log(1+f).   *   *	Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers - *		   n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last  - *		   20 bits (for VAX D format), or the last 21 bits ( for IEEE  + *		   n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last + *		   20 bits (for VAX D format), or the last 21 bits ( for IEEE   *		   double) is 0. This ensures n*ln2hi is exactly representable.   *		2. In step 1, f may not be representable. A correction term c   *	 	   for f is computed. It follows that the correction term for @@ -83,7 +83,7 @@ static char sccsid[] = "@(#)log1p.c	8.1 (Berkeley) 6/4/93";   *	only log1p(0)=0 is exact for finite argument.   *   * Accuracy: - *	log1p(x) returns the exact log(1+x) nearly rounded. In a test run  + *	log1p(x) returns the exact log(1+x) nearly rounded. In a test run   *	with 1,536,000 random arguments on a VAX, the maximum observed   *	error was .846 ulps (units in the last place).   * @@ -114,7 +114,7 @@ ic(sqrt2, 1.4142135623730951455E0,     0, 1.6A09E667F3BCD)  double log1p(x)  double x;  { -	const static double zero=0.0, negone= -1.0, one=1.0,  +	const static double zero=0.0, negone= -1.0, one=1.0,  		      half=1.0/2.0, small=1.0E-20;   /* 1+small == 1 */  	double z,s,t,c;  	int k; @@ -129,7 +129,7 @@ double x;  	   /* argument reduction */  	      if(copysign(x,one)<small) return(x);  	      k=logb(one+x); z=scalb(x,-k); t=scalb(one,-k); -	      if(z+t >= sqrt2 )  +	      if(z+t >= sqrt2 )  		  { k += 1 ; z *= half; t *= half; }  	      t += negone; x = z + t;  	      c = (t-x)+z ;		/* correction term for x */ @@ -162,9 +162,9 @@ double x;      /* end of if (finite(x)) */      /* log(-INF) is NaN */ -	else if(x<0)  +	else if(x<0)  	     return(zero/zero);      /* log(+INF) is INF */ -	else return(x);       +	else return(x);  } diff --git a/lib/libm/common_source/log__L.c b/lib/libm/common_source/log__L.c index c00158fa5172..207cb0d9e5d8 100644 --- a/lib/libm/common_source/log__L.c +++ b/lib/libm/common_source/log__L.c @@ -39,14 +39,14 @@ static char sccsid[] = "@(#)log__L.c	8.1 (Berkeley) 6/4/93";   *		LOG(1+X) - 2S			       X   * RETURN      ---------------  WHERE Z = S*S,  S = ------- , 0 <= Z <= .0294...   *		      S				     2 + X - *		      + *   * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)   * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS - * CODED IN C BY K.C. NG, 1/19/85;  + * CODED IN C BY K.C. NG, 1/19/85;   * REVISED BY K.C. Ng, 2/3/85, 4/16/85.   *   * Method : - *	1. Polynomial approximation: let s = x/(2+x).  + *	1. Polynomial approximation: let s = x/(2+x).   *	   Based on log(1+x) = log(1+s) - log(1-s)   *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,   * @@ -54,11 +54,11 @@ static char sccsid[] = "@(#)log__L.c	8.1 (Berkeley) 6/4/93";   *   *	       z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))   * - *	   where z=s*s. (See the listing below for Lk's values.) The  - *	   coefficients are obtained by a special Remez algorithm.  + *	   where z=s*s. (See the listing below for Lk's values.) The + *	   coefficients are obtained by a special Remez algorithm.   *   * Accuracy: - *	Assuming no rounding error, the maximum magnitude of the approximation  + *	Assuming no rounding error, the maximum magnitude of the approximation   *	error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)   *	for VAX D format.   * diff --git a/lib/libm/common_source/pow.c b/lib/libm/common_source/pow.c index cb4fc5029f4d..9e92985f8f97 100644 --- a/lib/libm/common_source/pow.c +++ b/lib/libm/common_source/pow.c @@ -35,17 +35,17 @@  static char sccsid[] = "@(#)pow.c	8.1 (Berkeley) 6/4/93";  #endif /* not lint */ -/* POW(X,Y)   - * RETURN X**Y  +/* POW(X,Y) + * RETURN X**Y   * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85;  + * CODED IN C BY K.C. NG, 1/8/85;   * REVISED BY K.C. NG on 7/10/85.   * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92.   * Required system supported functions: - *      scalb(x,n)       - *      logb(x)          - *	copysign(x,y)	 - *	finite(x)	 + *      scalb(x,n) + *      logb(x) + *	copysign(x,y) + *	finite(x)   *	drem(x,y)   *   * Required kernel functions: @@ -56,7 +56,7 @@ static char sccsid[] = "@(#)pow.c	8.1 (Berkeley) 6/4/93";   *	1. Compute and return log(x) in three pieces:   *		log(x) = n*ln2 + hi + lo,   *	   where n is an integer. - *	2. Perform y*log(x) by simulating muti-precision arithmetic and  + *	2. Perform y*log(x) by simulating muti-precision arithmetic and   *	   return the answer in three pieces:   *		y*log(x) = m*ln2 + hi + lo,   *	   where m is an integer. @@ -94,7 +94,7 @@ static char sccsid[] = "@(#)pow.c	8.1 (Berkeley) 6/4/93";   *			pow(integer,integer)   *	always returns the correct integer provided it is representable.   *	In a test run with 100,000 random arguments with 0 < x, y < 20.0 - *	on a VAX, the maximum observed error was 1.79 ulps (units in the  + *	on a VAX, the maximum observed error was 1.79 ulps (units in the   *	last place).   *   * Constants : @@ -123,7 +123,7 @@ const static double zero=0.0, one=1.0, two=2.0, negone= -1.0;  static double pow_P __P((double, double)); -double pow(x,y)  	 +double pow(x,y)  double x,y;  {  	double t; diff --git a/lib/libm/common_source/sinh.c b/lib/libm/common_source/sinh.c index 0516849cff13..075172c074e1 100644 --- a/lib/libm/common_source/sinh.c +++ b/lib/libm/common_source/sinh.c @@ -38,7 +38,7 @@ static char sccsid[] = "@(#)sinh.c	8.1 (Berkeley) 6/4/93";  /* SINH(X)   * RETURN THE HYPERBOLIC SINE OF X   * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85;  + * CODED IN C BY K.C. NG, 1/8/85;   * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85.   *   * Required system supported functions : @@ -50,14 +50,14 @@ static char sccsid[] = "@(#)sinh.c	8.1 (Berkeley) 6/4/93";   *   * Method :   *	1. reduce x to non-negative by sinh(-x) = - sinh(x). - *	2.  + *	2.   *   *	                                      expm1(x) + expm1(x)/(expm1(x)+1)   *	    0 <= x <= lnovfl     : sinh(x) := --------------------------------   *			       		                      2   *     lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow)   * lnovfl+ln2 <  x <  INF        :  overflow to INF - *	 + *   *   * Special cases:   *	sinh(x) is x if x is +INF, -INF, or NaN. @@ -112,7 +112,7 @@ double x;  	    {t=expm1(x); return(copysign((t+t/(one+t))*half,sign));}  	else if(x <= lnovfl+0.7) -		/* subtract x by ln(2^(max+1)) and return 2^max*exp(x)  +		/* subtract x by ln(2^(max+1)) and return 2^max*exp(x)  	    		to avoid unnecessary overflow */  	    return(copysign(scalb(one+expm1((x-mln2hi)-mln2lo),max),sign)); diff --git a/lib/libm/common_source/tanh.c b/lib/libm/common_source/tanh.c index d4923b3418bd..6813b55a3a70 100644 --- a/lib/libm/common_source/tanh.c +++ b/lib/libm/common_source/tanh.c @@ -38,7 +38,7 @@ static char sccsid[] = "@(#)tanh.c	8.1 (Berkeley) 6/4/93";  /* TANH(X)   * RETURN THE HYPERBOLIC TANGENT OF X   * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85;  + * CODED IN C BY K.C. NG, 1/8/85;   * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85.   *   * Required system supported functions : @@ -85,7 +85,7 @@ double x;  	sign=copysign(one,x);  	x=copysign(x,one); -	if(x < 22.0)  +	if(x < 22.0)  	    if( x > one )  		return(copysign(one-two/(expm1(x+x)+two),sign));  	    else if ( x > small ) | 
