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