- on non-ieee, rename gamma to tgamma, the 'true' gamma

- make gamma an alias of lgamma
- on ieee, add tgamma, based on gamma from non-ieee
- fixes for tgamma/lgamma/exp/log, esp. special cases (some from
freebsd);  properly raise invalid fp operations on vax
- also some general cleanup, ansification, man page (which was ok
jmc@)
- bump minor
this makes some ports using tgamma possible;  also consistifies
behavior across openbsd/ieee and openbsd/non-ieee, and other operating
systems
much thanks sthen@, johan@, steven@, Simon Kuhnle, Wiktor Izdebski
for testing
ok millert@

4 files changed, 123 insertions, 83 deletions

diff --git a/lib/libm/noieee_src/n_exp.c b/lib/libm/noieee_src/n_exp.c
index e1fef84bc1b..a51696b58d6 100644
--- a/lib/libm/noieee_src/n_exp.c
+++ b/lib/libm/noieee_src/n_exp.c
@@ -1,4 +1,4 @@
-/*      $NetBSD: n_exp.c,v 1.1 1995/10/10 23:36:44 ragge Exp $ */
+/*	$OpenBSD: n_exp.c,v 1.5 2008/06/11 20:53:27 martynas Exp $	*/
 /*
  * Copyright (c) 1985, 1993
  *	The Regents of the University of California.  All rights reserved.
@@ -35,21 +35,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:
- *	scalbn(x,n)	
- *	copysign(x,y)	
+ *	scalbn(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
@@ -111,10 +111,9 @@ ic(lnhuge, 7.1602103751842355450E2,    9,  1.6602B15B7ECF2)
 ic(lntiny,-7.5137154372698068983E2,    9, -1.77AF8EBEAE354)
 ic(invln2, 1.4426950408889633870E0,    0,  1.71547652B82FE)
 
-double exp(x)
-double x;
+double exp(double x)
 {
-	double  z,hi,lo,c;
+	double z, hi, lo, c;
 	int k;
 
 #if !defined(__vax__)&&!defined(tahoe)
@@ -140,7 +139,7 @@ double x;
 		}
 		/* end of x > lntiny */
 
-		else 
+		else
 		     /* exp(-big#) underflows to zero */
 		     if(finite(x))  return(scalbn(1.0,-5000));
 
@@ -149,17 +148,16 @@ double x;
 	}
 	/* end of x < lnhuge */
 
-	else 
+	else
 	/* exp(INF) is INF, exp(+big#) overflows to INF */
 	    return( finite(x) ?  scalbn(1.0,5000)  : x);
 }
 
 /* returns exp(r = x + c) for |c| < |x| with no overlap.  */
 
-double __exp__D(x, c)
-double x, c;
+double __exp__D(double x, double c)
 {
-	double  z,hi,lo, t;
+	double z, hi, lo;
 	int k;
 
 #if !defined(__vax__)&&!defined(tahoe)
@@ -185,7 +183,7 @@ double x, c;
 		}
 		/* end of x > lntiny */
 
-		else 
+		else
 		     /* exp(-big#) underflows to zero */
 		     if(finite(x))  return(scalbn(1.0,-5000));
 
@@ -194,7 +192,7 @@ double x, c;
 	}
 	/* end of x < lnhuge */
 
-	else 
+	else
 	/* exp(INF) is INF, exp(+big#) overflows to INF */
 	    return( finite(x) ?  scalbn(1.0,5000)  : x);
 }
diff --git a/lib/libm/noieee_src/n_lgamma.c b/lib/libm/noieee_src/n_lgamma.c
index 5eaed3d2f34..c874e8c6ba5 100644
--- a/lib/libm/noieee_src/n_lgamma.c
+++ b/lib/libm/noieee_src/n_lgamma.c
@@ -1,4 +1,4 @@
-/*      $NetBSD: n_lgamma.c,v 1.1 1995/10/10 23:36:56 ragge Exp $ */
+/*	$OpenBSD: n_lgamma.c,v 1.4 2008/06/11 20:53:27 martynas Exp $	*/
 /*-
  * Copyright (c) 1992, 1993
  *	The Regents of the University of California.  All rights reserved.
@@ -53,7 +53,7 @@ static char sccsid[] = "@(#)lgamma.c	8.2 (Berkeley) 11/30/93";
  *	x > 6:
  *		Use the asymptotic expansion (Stirling's Formula)
  *	0 < x < 6:
- *		Use gamma(x+1) = x*gamma(x) for argument reduction.
+ *		Use tgamma(x+1) = x*tgamma(x) for argument reduction.
  *		Use rational approximation in
  *		the range 1.2, 2.5
  *		Two approximations are used, one centered at the
@@ -66,12 +66,12 @@ static char sccsid[] = "@(#)lgamma.c	8.2 (Berkeley) 11/30/93";
  *	non-positive integer	returns +Inf.
  *	NaN			returns NaN
 */
-static int endian;
 #if defined(__vax__) || defined(tahoe)
 #define _IEEE		0
 /* double and float have same size exponent field */
 #define TRUNC(x)	x = (double) (float) (x)
 #else
+static int endian;
 #define _IEEE		1
 #define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
 #define infnan(x)	0.0
@@ -80,16 +80,16 @@ static int endian;
 static double small_lgam(double);
 static double large_lgam(double);
 static double neg_lgam(double);
-static double zero = 0.0, one = 1.0;
-int signgam;
+static const double one = 1.0;
+extern int signgam;
 
 #define UNDERFL (1e-1020 * 1e-1020)
 
 #define LEFT	(1.0 - (x0 + .25))
 #define RIGHT	(x0 - .218)
 /*
-/* Constants for approximation in [1.244,1.712]
-*/
+ * Constants for approximation in [1.244,1.712]
+ */
 #define x0	0.461632144968362356785
 #define x0_lo	-.000000000000000015522348162858676890521
 #define a0_hi	-0.12148629128932952880859
@@ -141,7 +141,9 @@ lgamma(double x)
 	double r;
 
 	signgam = 1;
+#if _IEEE
 	endian = ((*(int *) &one)) ? 1 : 0;
+#endif
 
 	if (!finite(x))
 		if (_IEEE)
@@ -161,11 +163,33 @@ lgamma(double x)
 		return (neg_lgam(x));
 }
 
+float
+lgammaf(float x)
+{
+	return lgamma(x);
+}
+
+/*
+ * The gamma() function performs identically to lgamma(), including
+ * the use of signgam.
+ */
+
+double
+gamma(double x)
+{
+	return lgamma(x);
+}
+
+float
+gammaf(float x)
+{
+	return lgammaf(x);
+}
+
 static double
 large_lgam(double x)
 {
 	double z, p, x1;
-	int i;
 	struct Double t, u, v;
 	u = __log__D(x);
 	u.a -= 1.0;
@@ -266,8 +290,7 @@ static double
 neg_lgam(double x)
 {
 	int xi;
-	double y, z, one = 1.0, zero = 0.0;
-	extern double gamma();
+	double y, z, zero = 0.0;
 
 	/* avoid destructive cancellation as much as possible */
 	if (x > -170) {
@@ -277,7 +300,7 @@ neg_lgam(double x)
 				return(one/zero);
 			else
 				return(infnan(ERANGE));
-		y = gamma(x);
+		y = tgamma(x);
 		if (y < 0)
 			y = -y, signgam = -1;
 		return (log(y));
diff --git a/lib/libm/noieee_src/n_log.c b/lib/libm/noieee_src/n_log.c
index ac367f72429..04f4a22bacc 100644
--- a/lib/libm/noieee_src/n_log.c
+++ b/lib/libm/noieee_src/n_log.c
@@ -1,4 +1,4 @@
-/*      $NetBSD: n_log.c,v 1.1 1995/10/10 23:36:57 ragge Exp $ */
+/*	$OpenBSD: n_log.c,v 1.4 2008/06/11 20:53:27 martynas Exp $	*/
 /*
  * Copyright (c) 1992, 1993
  *	The Regents of the University of California.  All rights reserved.
@@ -93,12 +93,12 @@ static char sccsid[] = "@(#)log.c	8.2 (Berkeley) 11/30/93";
  * Values for log(F) were generated using error < 10^-57 absolute
  * with the bc -l package.
 */
-static double	A1 = 	  .08333333333333178827;
-static double	A2 = 	  .01250000000377174923;
-static double	A3 =	 .002232139987919447809;
-static double	A4 =	.0004348877777076145742;
+static const double	A1 = 	  .08333333333333178827;
+static const double	A2 = 	  .01250000000377174923;
+static const double	A3 =	 .002232139987919447809;
+static const double	A4 =	.0004348877777076145742;
 
-static double logF_head[N+1] = {
+static const double logF_head[N+1] = {
 	0.,
 	.007782140442060381246,
 	.015504186535963526694,
@@ -230,7 +230,7 @@ static double logF_head[N+1] = {
 	.693147180560117703862
 };
 
-static double logF_tail[N+1] = {
+static const double logF_tail[N+1] = {
 	0.,
 	-.00000000000000543229938420049,
 	 .00000000000000172745674997061,
@@ -363,11 +363,7 @@ static double logF_tail[N+1] = {
 };
 
 double
-#ifdef _ANSI_SOURCE
 log(double x)
-#else
-log(x) double x;
-#endif
 {
 	int m, j;
 	double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0;
@@ -388,7 +384,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		*/
 
@@ -435,17 +431,13 @@ log(x) double x;
 
 /*
  * Extra precision variant, returning struct {double a, b;};
- * log(x) = a+b to 63 bits, with a is rounded to 26 bits.
+ * log(x) = a+b to 63 bits, with a rounded to 26 bits.
  */
 struct Double
-#ifdef _ANSI_SOURCE
 __log__D(double x)
-#else
-__log__D(x) double x;
-#endif
 {
 	int m, j;
-	double F, f, g, q, u, v, u2, one = 1.0;
+	double F, f, g, q, u, v, u2;
 	volatile double u1;
 	struct Double r;
 
diff --git a/lib/libm/noieee_src/n_gamma.c b/lib/libm/noieee_src/n_tgamma.c
index cce478274c7..de9ff444733 100644
--- a/lib/libm/noieee_src/n_gamma.c
+++ b/lib/libm/noieee_src/n_tgamma.c
@@ -1,4 +1,4 @@
-/*      $NetBSD: n_gamma.c,v 1.1 1995/10/10 23:36:50 ragge Exp $ */
+/*	$OpenBSD: n_tgamma.c,v 1.1 2008/06/11 20:53:27 martynas Exp $	*/
 /*-
  * Copyright (c) 1992, 1993
  *	The Regents of the University of California.  All rights reserved.
@@ -45,7 +45,7 @@ static char sccsid[] = "@(#)gamma.c	8.1 (Berkeley) 6/4/93";
 
 /* METHOD:
  * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x))
- * 	At negative integers, return +Inf, and set errno.
+ *	At negative integers, return +Inf and raise overflow.
  *
  * x < 6.5:
  *	Use argument reduction G(x+1) = xG(x) to reach the
@@ -62,11 +62,15 @@ static char sccsid[] = "@(#)gamma.c	8.1 (Berkeley) 6/4/93";
  *	avoid premature round-off.
  *
  * Special values:
- *	non-positive integer:	Set overflow trap; return +Inf;
- *	x > 171.63:		Set overflow trap; return +Inf;
- *	NaN: 			Set invalid trap;  return NaN
+ *	-Inf:			return 0 and raise invalid;
+ *	negative integer:	return +Inf and raise overflow;
+ *	other x ~< -177.79:	return 0 and raise underflow;
+ *	+-0:			return +Inf and raise overflow;
+ *	finite x ~> 171.63:	return +Inf and raise div-by-0;
+ *	+Inf:			return +Inf and raise div-by-0;
+ *	NaN:			return 0 and raise invalid.
  *
- * Accuracy: Gamma(x) is accurate to within
+ * Accuracy: tgamma(x) is accurate to within
  *	x > 0:  error provably < 0.9ulp.
  *	Maximum observed in 1,000,000 trials was .87ulp.
  *	x < 0:
@@ -118,7 +122,7 @@ static struct Double ratfun_gam(double, double);
 #define Pa7	-1.44705562421428915453880392761e-02
 
 static const double zero = 0., one = 1.0, tiny = 1e-300;
-static int endian;
+
 /*
  * TRUNC sets trailing bits in a floating-point number to zero.
  * is a temporary variable.
@@ -127,21 +131,26 @@ static int endian;
 #define _IEEE		0
 #define TRUNC(x)	x = (double) (float) (x)
 #else
+static int endian;
 #define _IEEE		1
 #define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
 #define infnan(x)	0.0
 #endif
 
 double
-gamma(x)
-	double x;
+tgamma(double x)
 {
 	struct Double u;
+#if _IEEE
 	endian = (*(int *) &one) ? 1 : 0;
+#endif
 
 	if (x >= 6) {
 		if(x > 171.63)
-			return(one/zero);
+			if (_IEEE)
+				return (x/zero);
+			else
+				return (infnan(ERANGE));
 		u = large_gam(x);
 		return(__exp__D(u.a, u.b));
 	} else if (x >= 1.0 + LEFT + x0)
@@ -149,28 +158,43 @@ gamma(x)
 	else if (x > 1.e-17)
 		return (smaller_gam(x));
 	else if (x > -1.e-17) {
-		if (x == 0.0)
-			if (!_IEEE) return (infnan(ERANGE));
-			else return (one/x);
-		one+1e-20;		/* Raise inexact flag. */
+		if (x == 0.0) {
+			if (!_IEEE)
+				return (infnan(ERANGE));
+		} else {
+			u.a = one - tiny;	/* raise inexact */
+		}
 		return (one/x);
 	} else if (!finite(x)) {
 		if (_IEEE)		/* x = NaN, -Inf */
-			return (x*x);
+			return (x - x);
 		else
 			return (infnan(EDOM));
 	 } else
 		return (neg_gam(x));
 }
+
+/*
+ * We simply call tgamma() rather than bloating the math library
+ * with a float-optimized version of it.  The reason is that tgammaf()
+ * is essentially useless, since the function is superexponential
+ * and floats have very limited range.  -- das@freebsd.org
+ */
+
+float
+tgammaf(float x)
+{
+	return tgamma(x);
+}
+
 /*
  * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
  */
+
 static struct Double
-large_gam(x)
-	double x;
+large_gam(double x)
 {
 	double z, p;
-	int i;
 	struct Double t, u, v;
 
 	z = one/(x*x);
@@ -191,15 +215,16 @@ large_gam(x)
 	u.b += lns2pi_hi; u.b += t.b;
 	return (u);
 }
+
 /*
  * Good to < 1 ulp.  (provably .90 ulp; .87 ulp on 1,000,000 runs.)
  * It also has correct monotonicity.
  */
+
 static double
-small_gam(x)
-	double x;
+small_gam(double x)
 {
-	double y, ym1, t, x1;
+	double y, ym1, t;
 	struct Double yy, r;
 	y = x - one;
 	ym1 = y - one;
@@ -220,18 +245,19 @@ small_gam(x)
 		TRUNC(r.a);
 		r.b += (t - r.a);
 	}
-	/* Return r*gamma(y). */
+	/* Return r*tgamma(y). */
 	yy = ratfun_gam(y - x0, 0);
 	y = r.b*(yy.a + yy.b) + r.a*yy.b;
 	y += yy.a*r.a;
 	return (y);
 }
+
 /*
  * Good on (0, 1+x0+LEFT].  Accurate to 1ulp.
  */
+
 static double
-smaller_gam(x)
-	double x;
+smaller_gam(double x)
 {
 	double t, d;
 	struct Double r, xx;
@@ -255,14 +281,14 @@ smaller_gam(x)
 	r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b;
 	return (d + r.a/x);
 }
+
 /*
  * returns (z+c)^2 * P(z)/Q(z) + a0
  */
+
 static struct Double
-ratfun_gam(z, c)
-	double z, c;
+ratfun_gam(double z, double c)
 {
-	int i;
 	double p, q;
 	struct Double r, t;
 
@@ -287,21 +313,22 @@ ratfun_gam(z, c)
 }
 
 static double
-neg_gam(x)
-	double x;
+neg_gam(double x)
 {
 	int sgn = 1;
 	struct Double lg, lsine;
 	double y, z;
 
-	y = floor(x + .5);
+	y = ceil(x);
 	if (y == x)		/* Negative integer. */
-		if(!_IEEE)
-			return (infnan(ERANGE));
+		if (_IEEE)
+			return ((x - x) / zero);
 		else
-			return (one/zero);
-	z = fabs(x - y);
-	y = .5*ceil(x);
+			return (infnan(ERANGE));
+	z = y - x;
+	if (z > 0.5)
+		z = one - z;
+	y = 0.5 * y;
 	if (y == ceil(y))
 		sgn = -1;
 	if (z < .25)
@@ -325,9 +352,9 @@ neg_gam(x)
 	}
 	y = one-x;
 	if (one-y == x)
-		y = gamma(y);
+		y = tgamma(y);
 	else		/* 1-x is inexact */
-		y = -x*gamma(-x);
+		y = -x*tgamma(-x);
 	if (sgn < 0) y = -y;
 	return (M_PI / (y*z));
 }