diff options
27 files changed, 270 insertions, 243 deletions
diff --git a/lib/libm/noieee_src/mathimpl.h b/lib/libm/noieee_src/mathimpl.h index 7fd4fa8a930..0439bf710b6 100644 --- a/lib/libm/noieee_src/mathimpl.h +++ b/lib/libm/noieee_src/mathimpl.h @@ -1,3 +1,4 @@ +/* $OpenBSD: mathimpl.h,v 1.7 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: mathimpl.h,v 1.1 1995/10/10 23:36:31 ragge Exp $ */ /* * Copyright (c) 1988, 1993 diff --git a/lib/libm/noieee_src/n_acosh.c b/lib/libm/noieee_src/n_acosh.c index 652624c3248..07faacead2f 100644 --- a/lib/libm/noieee_src/n_acosh.c +++ b/lib/libm/noieee_src/n_acosh.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_acosh.c,v 1.6 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_acosh.c,v 1.1 1995/10/10 23:36:33 ragge Exp $ */ /* * Copyright (c) 1985, 1993 @@ -45,10 +46,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. * @@ -57,7 +58,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 . @@ -86,14 +87,14 @@ ic(ln2lo, 1.9082149292705877000E-10,-33, 1.A39EF35793C76) double acosh(x) double x; -{ +{ double t,big=1.E20; /* big+1==big */ if (isnan(x)) return (x); /* 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/noieee_src/n_asincos.c b/lib/libm/noieee_src/n_asincos.c index 7e99c1c66f6..b8d6567f9cd 100644 --- a/lib/libm/noieee_src/n_asincos.c +++ b/lib/libm/noieee_src/n_asincos.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_asincos.c,v 1.6 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_asincos.c,v 1.1 1995/10/10 23:36:34 ragge Exp $ */ /* * Copyright (c) 1985, 1993 @@ -42,12 +43,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: @@ -56,22 +57,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))); * @@ -79,7 +80,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. */ @@ -99,7 +100,7 @@ asin(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))); } } @@ -114,9 +115,9 @@ asin(x) * sqrt(x) * * Required kernel function: - * atan2(y,x) + * atan2(y,x) * - * Method : + * Method: * ________ * / 1 - x * acos(x) = 2*atan2( / -------- , 1 ) . @@ -128,22 +129,22 @@ asin(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))); * @@ -151,7 +152,7 @@ asin(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/noieee_src/n_asinh.c b/lib/libm/noieee_src/n_asinh.c index c5e7f867922..55e8bd9b53e 100644 --- a/lib/libm/noieee_src/n_asinh.c +++ b/lib/libm/noieee_src/n_asinh.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_asinh.c,v 1.6 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_asinh.c,v 1.1 1995/10/10 23:36:35 ragge Exp $ */ /* * Copyright (c) 1985, 1993 @@ -46,16 +47,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: @@ -81,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 (isnan(x)) return (x); - 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/noieee_src/n_atan.c b/lib/libm/noieee_src/n_atan.c index 0d3d8ff40e5..c0af9a97122 100644 --- a/lib/libm/noieee_src/n_atan.c +++ b/lib/libm/noieee_src/n_atan.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_atan.c,v 1.4 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_atan.c,v 1.1 1995/10/10 23:36:36 ragge Exp $ */ /* * Copyright (c) 1985, 1993 @@ -38,32 +39,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))). * @@ -71,7 +72,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/noieee_src/n_atan2.c b/lib/libm/noieee_src/n_atan2.c index e8f453f29d3..3522f0865d3 100644 --- a/lib/libm/noieee_src/n_atan2.c +++ b/lib/libm/noieee_src/n_atan2.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_atan2.c,v 1.1 1995/10/10 23:36:37 ragge Exp $ */ +/* $OpenBSD: n_atan2.c,v 1.7 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_atan2.c,v 1.1 1995/10/10 23:36:37 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -35,21 +36,21 @@ static char sccsid[] = "@(#)atan2.c 8.1 (Berkeley) 6/4/93"; /* ATAN2(Y,X) * RETURN ARG (X+iY) * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; + * CODED IN C BY K.C. NG, 1/8/85; * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85. * * Required system supported functions : * copysign(x,y) * scalbn(x,y) * logb(x) - * + * * Method : * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). - * 2. Reduce x to positive by (if x and y are unexceptional): + * 2. Reduce x to positive by (if x and y are unexceptional): * ARG (x+iy) = arctan(y/x) ... if x > 0, * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, - * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument - * is further reduced to one of the following intervals and the + * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument + * is further reduced to one of the following intervals and the * arctangent of y/x is evaluated by the corresponding formula: * * [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) @@ -73,32 +74,32 @@ static char sccsid[] = "@(#)atan2.c 8.1 (Berkeley) 6/4/93"; * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; * * Accuracy: - * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, + * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, * where * * in decimal: - * pi = 3.141592653589793 23846264338327 ..... + * pi = 3.141592653589793 23846264338327 ..... * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , + * 56 bits PI = 3.141592653589793 227020265 ..... , * * in hexadecimal: * pi = 3.243F6A8885A308D313198A2E.... * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * + * * In a test run with 356,000 random argument on [-1,1] * [-1,1] on a * VAX, the maximum observed error was 1.41 ulps (units of the last place) * compared with (PI/pi)*(the exact ARG(x+iy)). * * Note: * We use machine PI (the true pi rounded) in place of the actual - * value of pi for all the trig and inverse trig functions. In general, - * if trig is one of sin, cos, tan, then computed trig(y) returns the - * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig - * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the + * value of pi for all the trig and inverse trig functions. In general, + * if trig is one of sin, cos, tan, then computed trig(y) returns the + * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig + * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the * trig functions have period PI, and trig(arctrig(x)) returns x for * all critical values x. - * + * * Constants: * The hexadecimal values are the intended ones for the following constants. * The decimal values may be used, provided that the compiler will convert @@ -172,7 +173,7 @@ ic(a11, 1.6438029044759730479E-2 , -6, 1.0D52174A1BB54) double atan2(y,x) double y,x; -{ +{ static const double zero=0, one=1, small=1.0E-9, big=1.0E18; double t,z,signy,signx,hi,lo; int k,m; @@ -184,8 +185,8 @@ double y,x; return (y); /* copy down the sign of y and x */ - signy = copysign(one,y) ; - signx = copysign(one,x) ; + signy = copysign(one,y) ; + signx = copysign(one,x) ; /* if x is 1.0, goto begin */ if(x==1) { y=copysign(y,one); t=y; if(finite(t)) goto begin;} @@ -195,10 +196,10 @@ double y,x; /* when x = 0 */ if(x==zero) return(copysign(PIo2,signy)); - + /* when x is INF */ if(!finite(x)) - if(!finite(y)) + if(!finite(y)) return(copysign((signx==one)?PIo4:3*PIo4,signy)); else return(copysign((signx==one)?zero:PI,signy)); @@ -207,43 +208,43 @@ double y,x; if(!finite(y)) return(copysign(PIo2,signy)); /* compute y/x */ - x=copysign(x,one); - y=copysign(y,one); - if((m=(k=logb(y))-logb(x)) > 60) t=big+big; + x=copysign(x,one); + y=copysign(y,one); + if((m=(k=logb(y))-logb(x)) > 60) t=big+big; else if(m < -80 ) t=y/x; else { t = y/x ; y = scalbn(y,-k); x=scalbn(x,-k); } /* begin argument reduction */ begin: - if (t < 2.4375) { + if (t < 2.4375) { /* truncate 4(t+1/16) to integer for branching */ k = 4 * (t+0.0625); switch (k) { /* t is in [0,7/16] */ - case 0: + case 0: case 1: - if (t < small) + if (t < small) { big + small ; /* raise inexact flag */ return (copysign((signx>zero)?t:PI-t,signy)); } hi = zero; lo = zero; break; /* t is in [7/16,11/16] */ - case 2: + case 2: hi = athfhi; lo = athflo; z = x+x; t = ( (y+y) - x ) / ( z + y ); break; /* t is in [11/16,19/16] */ - case 3: + case 3: case 4: hi = PIo4; lo = zero; t = ( y - x ) / ( x + y ); break; /* t is in [19/16,39/16] */ - default: + default: hi = at1fhi; lo = at1flo; z = y-x; y=y+y+y; t = x+x; t = ( (z+z)-x ) / ( t + y ); break; @@ -251,7 +252,7 @@ begin: } /* end of if (t < 2.4375) */ - else + else { hi = PIo2; lo = zero; @@ -259,7 +260,7 @@ begin: if (t <= big) t = - x / y; /* t is in [big, INF] */ - else + else { big+small; /* raise inexact flag */ t = zero; } } diff --git a/lib/libm/noieee_src/n_atanh.c b/lib/libm/noieee_src/n_atanh.c index f5dd6e327d5..ad4dc4a668d 100644 --- a/lib/libm/noieee_src/n_atanh.c +++ b/lib/libm/noieee_src/n_atanh.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_atanh.c,v 1.1 1995/10/10 23:36:38 ragge Exp $ */ +/* $OpenBSD: n_atanh.c,v 1.5 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_atanh.c,v 1.1 1995/10/10 23:36:38 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -35,14 +36,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 + * Method: + * Return * 1 2x x * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) * 2 1 - x 1 - x diff --git a/lib/libm/noieee_src/n_cabs.c b/lib/libm/noieee_src/n_cabs.c index 695306117f8..5c15a847211 100644 --- a/lib/libm/noieee_src/n_cabs.c +++ b/lib/libm/noieee_src/n_cabs.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_cabs.c,v 1.1 1995/10/10 23:36:39 ragge Exp $ */ +/* $OpenBSD: n_cabs.c,v 1.7 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_cabs.c,v 1.1 1995/10/10 23:36:39 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -35,7 +36,7 @@ static char sccsid[] = "@(#)cabs.c 8.1 (Berkeley) 6/4/93"; /* HYPOT(X,Y) * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 11/28/84; + * CODED IN C BY K.C. NG, 11/28/84; * REVISED BY K.C. NG, 7/12/85. * * Required system supported functions : @@ -49,16 +50,16 @@ static char sccsid[] = "@(#)cabs.c 8.1 (Berkeley) 6/4/93"; * y if y > x (hence x is never smaller than y). * 2. Hypot(x,y) is computed by: * Case I, x/y > 2 - * + * * y * hypot = x + ----------------------------- * 2 * sqrt ( 1 + [x/y] ) + x/y * - * Case II, x/y <= 2 + * Case II, x/y <= 2 * y * hypot = x + -------------------------------------------------- - * 2 + * 2 * [x/y] - 2 * (sqrt(2)+1) + (x-y)/y + ----------------------------- * 2 @@ -106,7 +107,7 @@ double hypot(x,y) double x, y; { - static const double zero=0, one=1, + static const double zero=0, one=1, small=1.0E-18; /* fl(1+small)==1 */ static const ibig=30; /* fl(1+2**(2*ibig))==1 */ double t,r; @@ -114,15 +115,15 @@ double x, y; if(finite(x)) if(finite(y)) - { + { x=copysign(x,one); y=copysign(y,one); - if(y > x) + if(y > x) { t=x; x=y; y=t; } if(x == zero) return(zero); if(y == zero) return(x); exp= logb(x); - if(exp-(int)logb(y) > ibig ) + if(exp-(int)logb(y) > ibig ) /* raise inexact flag and return |x| */ { one+small; return(x); } @@ -143,7 +144,7 @@ double x, y; else if(y==y) /* y is +-INF */ return(copysign(y,one)); - else + else return(y); /* y is NaN and x is finite */ else if(x==x) /* x is +-INF */ @@ -197,16 +198,16 @@ double x, y; if(finite(x)) if(finite(y)) - { + { x=copysign(x,one); y=copysign(y,one); - if(y > x) + if(y > x) { temp=x; x=y; y=temp; } if(x == zero) return(zero); if(y == zero) return(x); exp= logb(x); x=scalbn(x,-exp); - if(exp-(int)logb(y) > ibig ) + if(exp-(int)logb(y) > ibig ) /* raise inexact flag and return |x| */ { one+small; return(scalbn(x,exp)); } else y=scalbn(y,-exp); @@ -215,7 +216,7 @@ double x, y; else if(y==y) /* y is +-INF */ return(copysign(y,one)); - else + else return(y); /* y is NaN and x is finite */ else if(x==x) /* x is +-INF */ diff --git a/lib/libm/noieee_src/n_cbrt.c b/lib/libm/noieee_src/n_cbrt.c index f0aafd54f38..03f23d4e657 100644 --- a/lib/libm/noieee_src/n_cbrt.c +++ b/lib/libm/noieee_src/n_cbrt.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_cbrt.c,v 1.4 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_cbrt.c,v 1.1 1995/10/10 23:36:40 ragge Exp $ */ /* * Copyright (c) 1985, 1993 @@ -47,8 +48,8 @@ static char sccsid[] = "@(#)cbrt.c 8.1 (Berkeley) 6/4/93"; * long interger at the address of a floating point number will be the * leading 32 bits of that floating point number (i.e., sign, exponent, * and the 20 most significant bits). - * On a National machine, it has different ordering; therefore, this code - * must be compiled with flag -DNATIONAL. + * On a National machine, it has different ordering; therefore, this code + * must be compiled with flag -DNATIONAL. */ #if !defined(__vax__)&&!defined(tahoe) @@ -62,7 +63,7 @@ static const double F= 45./28., G= 5./14.; -double cbrt(x) +double cbrt(x) double x; { double r,s,t=0.0,w; @@ -90,15 +91,15 @@ double x; t*=x; pt[n0]=pt[n0]/3+B2; } else - pt[n0]=px[n0]/3+B1; + pt[n0]=px[n0]/3+B1; /* new cbrt to 23 bits, may be implemented in single precision */ r=t*t/x; s=C+r*t; - t*=G+F/(s+E+D/s); + t*=G+F/(s+E+D/s); - /* chopped to 20 bits and make it larger than cbrt(x) */ + /* chopped to 20 bits and make it larger than cbrt(x) */ pt[n1]=0; pt[n0]+=0x00000001; diff --git a/lib/libm/noieee_src/n_cosh.c b/lib/libm/noieee_src/n_cosh.c index cd028e3fa88..a95dd73d88a 100644 --- a/lib/libm/noieee_src/n_cosh.c +++ b/lib/libm/noieee_src/n_cosh.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_cosh.c,v 1.1 1995/10/10 23:36:42 ragge Exp $ */ +/* $OpenBSD: n_cosh.c,v 1.7 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_cosh.c,v 1.1 1995/10/10 23:36:42 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -35,7 +36,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 : @@ -43,20 +44,20 @@ static char sccsid[] = "@(#)cosh.c 8.1 (Berkeley) 6/4/93"; * scalbn(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 @@ -104,7 +105,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; @@ -113,19 +114,19 @@ double x; return (x); 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(scalbn(exp((x-mln2hi)-mln2lo), max)); + return(scalbn(exp((x-mln2hi)-mln2lo), max)); - else + else return(exp(x)*half); /* for large x, cosh(x)=exp(x)/2 */ } diff --git a/lib/libm/noieee_src/n_erf.c b/lib/libm/noieee_src/n_erf.c index 38d301793a1..3281985a8d9 100644 --- a/lib/libm/noieee_src/n_erf.c +++ b/lib/libm/noieee_src/n_erf.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_erf.c,v 1.5 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_erf.c,v 1.1 1995/10/10 23:36:43 ragge Exp $ */ /*- * Copyright (c) 1992, 1993 @@ -68,13 +69,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 +97,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 +106,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 +140,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 +182,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 +320,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; @@ -352,7 +353,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/noieee_src/n_exp__E.c b/lib/libm/noieee_src/n_exp__E.c index b55d8e76935..b00d0b8e562 100644 --- a/lib/libm/noieee_src/n_exp__E.c +++ b/lib/libm/noieee_src/n_exp__E.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_exp__E.c,v 1.1 1995/10/10 23:36:45 ragge Exp $ */ +/* $OpenBSD: n_exp__E.c,v 1.5 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_exp__E.c,v 1.1 1995/10/10 23:36:45 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -38,7 +39,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) @@ -47,12 +48,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 @@ -73,7 +74,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: @@ -118,7 +119,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/noieee_src/n_expm1.c b/lib/libm/noieee_src/n_expm1.c index fc2ebe51a3e..8dca54117b5 100644 --- a/lib/libm/noieee_src/n_expm1.c +++ b/lib/libm/noieee_src/n_expm1.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_expm1.c,v 1.1 1995/10/10 23:36:46 ragge Exp $ */ +/* $OpenBSD: n_expm1.c,v 1.7 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_expm1.c,v 1.1 1995/10/10 23:36:46 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -35,36 +36,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: - * scalbn(x,n) - * copysign(x,y) + * scalbn(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; @@ -106,7 +107,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) @@ -123,13 +124,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);} @@ -140,17 +141,17 @@ double x; { x=one-scalbn(one,-k); z += __exp__E(z,c);} else if(k<100) { x = __exp__E(z,c)-scalbn(one,-k); x+=z; z=one;} - else + else { x = __exp__E(z,c)+z; z=one;} - return (scalbn(x+z,k)); + return (scalbn(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 */ @@ -158,7 +159,7 @@ double x; } /* end of x < lnhuge */ - else + else /* expm1(INF) is INF, expm1(+big#) overflows to INF */ return( finite(x) ? scalbn(one,5000) : x); } diff --git a/lib/libm/noieee_src/n_floor.c b/lib/libm/noieee_src/n_floor.c index 05b29359e19..731263e53aa 100644 --- a/lib/libm/noieee_src/n_floor.c +++ b/lib/libm/noieee_src/n_floor.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_floor.c,v 1.1 1995/10/10 23:36:48 ragge Exp $ */ +/* $OpenBSD: n_floor.c,v 1.6 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_floor.c,v 1.1 1995/10/10 23:36:48 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. diff --git a/lib/libm/noieee_src/n_fmod.c b/lib/libm/noieee_src/n_fmod.c index b4b1c0f8860..c7e5eb2b336 100644 --- a/lib/libm/noieee_src/n_fmod.c +++ b/lib/libm/noieee_src/n_fmod.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_fmod.c,v 1.5 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_fmod.c,v 1.1 1995/10/10 23:36:49 ragge Exp $ */ /* * Copyright (c) 1989, 1993 diff --git a/lib/libm/noieee_src/n_j0.c b/lib/libm/noieee_src/n_j0.c index c916cb93086..33b8e8f5040 100644 --- a/lib/libm/noieee_src/n_j0.c +++ b/lib/libm/noieee_src/n_j0.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_j0.c,v 1.5 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_j0.c,v 1.1 1995/10/10 23:36:52 ragge Exp $ */ /*- * Copyright (c) 1992, 1993 @@ -43,18 +44,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 * @@ -81,20 +82,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. @@ -120,7 +121,7 @@ static char sccsid[] = "@(#)j0.c 8.2 (Berkeley) 11/30/93"; static double pzero(double), qzero(double); -static double +static double huge = 1e300, zero = 0.0, one = 1.0, @@ -137,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; @@ -200,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; @@ -344,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/noieee_src/n_j1.c b/lib/libm/noieee_src/n_j1.c index 3e251104acb..f399afe306b 100644 --- a/lib/libm/noieee_src/n_j1.c +++ b/lib/libm/noieee_src/n_j1.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_j1.c,v 1.4 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_j1.c,v 1.1 1995/10/10 23:36:53 ragge Exp $ */ /*- * Copyright (c) 1992, 1993 @@ -43,18 +44,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 * @@ -82,16 +83,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, @@ -120,7 +121,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, @@ -140,7 +141,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; @@ -203,7 +204,7 @@ static double v0[5] = { 1.665592462079920695971450872592458916421e-0011, }; -double y1(x) +double y1(x) double x; { double z, s, c, ss, cc, u, v; @@ -252,10 +253,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])))); @@ -349,7 +350,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/noieee_src/n_jn.c b/lib/libm/noieee_src/n_jn.c index 402aadda382..b1ba14cc364 100644 --- a/lib/libm/noieee_src/n_jn.c +++ b/lib/libm/noieee_src/n_jn.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_jn.c,v 1.4 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_jn.c,v 1.1 1995/10/10 23:36:54 ragge Exp $ */ /*- * Copyright (c) 1992, 1993 @@ -43,18 +44,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 * @@ -66,7 +67,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. @@ -85,7 +86,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> @@ -117,7 +118,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; } @@ -131,10 +132,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 @@ -151,7 +152,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++){ @@ -162,7 +163,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 */ @@ -177,14 +178,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 @@ -200,9 +201,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; @@ -251,7 +252,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; @@ -272,10 +273,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/noieee_src/n_log10.c b/lib/libm/noieee_src/n_log10.c index 4b6d5eb2578..111704542c7 100644 --- a/lib/libm/noieee_src/n_log10.c +++ b/lib/libm/noieee_src/n_log10.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_log10.c,v 1.1 1995/10/10 23:36:58 ragge Exp $ */ +/* $OpenBSD: n_log10.c,v 1.5 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_log10.c,v 1.1 1995/10/10 23:36:58 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -35,9 +36,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) * @@ -49,12 +50,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/noieee_src/n_log1p.c b/lib/libm/noieee_src/n_log1p.c index 1f226ad23db..dd60a2dd5ac 100644 --- a/lib/libm/noieee_src/n_log1p.c +++ b/lib/libm/noieee_src/n_log1p.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_log1p.c,v 1.1 1995/10/10 23:37:00 ragge Exp $ */ +/* $OpenBSD: n_log1p.c,v 1.7 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_log1p.c,v 1.1 1995/10/10 23:37:00 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -32,24 +33,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: - * scalbn(x,n) + * scalbn(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) @@ -62,11 +63,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 @@ -80,7 +81,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). * @@ -112,7 +113,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; @@ -126,7 +127,7 @@ double x; /* argument reduction */ if(copysign(x,one)<small) return(x); k=logb(one+x); z=scalbn(x,-k); t=scalbn(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 */ @@ -159,9 +160,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/noieee_src/n_log__L.c b/lib/libm/noieee_src/n_log__L.c index 895760a405e..02b13b1e10e 100644 --- a/lib/libm/noieee_src/n_log__L.c +++ b/lib/libm/noieee_src/n_log__L.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_log__L.c,v 1.1 1995/10/10 23:37:01 ragge Exp $ */ +/* $OpenBSD: n_log__L.c,v 1.5 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_log__L.c,v 1.1 1995/10/10 23:37:01 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -36,14 +37,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 + ....., * @@ -51,11 +52,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/noieee_src/n_pow.c b/lib/libm/noieee_src/n_pow.c index 3812a0c03cc..3043740e511 100644 --- a/lib/libm/noieee_src/n_pow.c +++ b/lib/libm/noieee_src/n_pow.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_pow.c,v 1.1 1995/10/10 23:37:02 ragge Exp $ */ +/* $OpenBSD: n_pow.c,v 1.8 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_pow.c,v 1.1 1995/10/10 23:37:02 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -32,17 +33,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: - * scalbn(x,n) - * logb(x) - * copysign(x,y) - * finite(x) + * scalbn(x,n) + * logb(x) + * copysign(x,y) + * finite(x) * drem(x,y) * * Required kernel functions: @@ -53,7 +54,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. @@ -91,7 +92,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 : @@ -120,7 +121,7 @@ const static double zero=0.0, one=1.0, two=2.0, negone= -1.0; static double pow_P(double, double); -double pow(x,y) +double pow(x,y) double x,y; { double t; diff --git a/lib/libm/noieee_src/n_sincos.c b/lib/libm/noieee_src/n_sincos.c index 6b8796439bc..6bf91bef357 100644 --- a/lib/libm/noieee_src/n_sincos.c +++ b/lib/libm/noieee_src/n_sincos.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_sincos.c,v 1.4 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_sincos.c,v 1.1 1995/10/10 23:37:04 ragge Exp $ */ /* * Copyright (c) 1987, 1993 @@ -66,7 +67,7 @@ double x; } double -cos(x) +cos(x) double x; { double a,c,z,s = 1.0; @@ -82,7 +83,7 @@ double x; } else { /* ... in [PI/4,3PI/4] */ a = PIo2-a; - return a+a*sin__S(a*a); /* rtn. S(PI/2-|x|) */ + return a+a*sin__S(a*a); /* rtn. S(PI/2-|x|) */ } } if (a < small) { diff --git a/lib/libm/noieee_src/n_sinh.c b/lib/libm/noieee_src/n_sinh.c index 77dd4d921a3..3747477922e 100644 --- a/lib/libm/noieee_src/n_sinh.c +++ b/lib/libm/noieee_src/n_sinh.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_sinh.c,v 1.1 1995/10/10 23:37:05 ragge Exp $ */ +/* $OpenBSD: n_sinh.c,v 1.7 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_sinh.c,v 1.1 1995/10/10 23:37:05 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. @@ -35,7 +36,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 : @@ -47,14 +48,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. @@ -111,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(scalbn(one+expm1((x-mln2hi)-mln2lo),max),sign)); diff --git a/lib/libm/noieee_src/n_support.c b/lib/libm/noieee_src/n_support.c index fa5d186112e..4b8549e18ea 100644 --- a/lib/libm/noieee_src/n_support.c +++ b/lib/libm/noieee_src/n_support.c @@ -1,4 +1,5 @@ -/* $NetBSD: n_support.c,v 1.1 1995/10/10 23:37:06 ragge Exp $ */ +/* $OpenBSD: n_support.c,v 1.8 2008/06/12 22:43:36 martynas Exp $ */ +/* $NetBSD: n_support.c,v 1.1 1995/10/10 23:37:06 ragge Exp $ */ /* * Copyright (c) 1985, 1993 * The Regents of the University of California. All rights reserved. diff --git a/lib/libm/noieee_src/n_tan.c b/lib/libm/noieee_src/n_tan.c index f39a1958397..edc483636cd 100644 --- a/lib/libm/noieee_src/n_tan.c +++ b/lib/libm/noieee_src/n_tan.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_tan.c,v 1.4 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_tan.c,v 1.1 1995/10/10 23:37:07 ragge Exp $ */ /* * Copyright (c) 1987, 1993 @@ -36,7 +37,7 @@ static char sccsid[] = "@(#)tan.c 8.1 (Berkeley) 6/4/93"; #include "mathimpl.h" double -tan(x) +tan(x) double x; { double a,z,ss,cc,c; diff --git a/lib/libm/noieee_src/n_tanh.c b/lib/libm/noieee_src/n_tanh.c index bf41f6c3b90..ab2015da174 100644 --- a/lib/libm/noieee_src/n_tanh.c +++ b/lib/libm/noieee_src/n_tanh.c @@ -1,3 +1,4 @@ +/* $OpenBSD: n_tanh.c,v 1.6 2008/06/12 22:43:36 martynas Exp $ */ /* $NetBSD: n_tanh.c,v 1.1 1995/10/10 23:37:08 ragge Exp $ */ /* * Copyright (c) 1985, 1993 @@ -37,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 : @@ -82,7 +83,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 ) |