libm: Fix 'gamma' and 'gammaf' functions. Clean up other gamma code. [v2]

The current gamma, gamma_r, gammaf and gammaf_r functions return
|gamma(x)| instead of ln(|gamma(x)|) due to a change made back in 2002
to the __ieee754_gamma_r implementation. This patch fixes that, making
all of these functions map too their lgamma equivalents.

To fix the underlying bug, the __ieee754_gamma functions have been
changed to return gamma(x), removing the _r variants as those are no
longer necessary. Their names have been changed to __ieee754_tgamma to
avoid potential confusion from users.

Now that the __ieee754_tgamma functions return the correctly signed
value, the tgamma functions have been modified to use them.

libm.a now exposes the following gamma functions:

    ln(|gamma(x)|):

	__ieee754_lgamma_r
	__ieee754_lgammaf_r

	lgamma
	lgamma_r
	gamma
	gamma_r

	lgammaf
	lgammaf_r
	gammaf
	gammaf_r

	lgammal	(on machines where long double is double)

    gamma(x):

	__ieee754_tgamma
	__ieee754_tgammaf
	tgamma
	tgammaf
	tgammal (on machines where long double is double)

Additional aliases for any of the above functions can be added if
necessary; in particular, I'm not sure if we need to include
__ieee754_gamma*_r functions (which would return ln(|(gamma(x)|).

Signed-off-by: Keith Packard <keithp@keithp.com>

----

v2:
	Switch commit message to ASCII

12 files changed, 53 insertions, 109 deletions

diff --git a/newlib/libm/math/Makefile.am b/newlib/libm/math/Makefile.am
index e745159ae..05cf87a4d 100644
--- a/newlib/libm/math/Makefile.am
+++ b/newlib/libm/math/Makefile.am
@@ -8,7 +8,7 @@ src = 	k_standard.c k_rem_pio2.c \
 	k_cos.c k_sin.c k_tan.c \
 	e_acos.c e_acosh.c e_asin.c e_atan2.c \
 	e_atanh.c e_cosh.c e_exp.c e_fmod.c \
-	er_gamma.c e_hypot.c e_j0.c \
+	e_tgamma.c e_hypot.c e_j0.c \
 	e_j1.c e_jn.c er_lgamma.c \
 	e_log.c e_log10.c e_pow.c e_rem_pio2.c e_remainder.c \
 	e_scalb.c e_sinh.c e_sqrt.c \
@@ -31,7 +31,7 @@ fsrc =	kf_rem_pio2.c \
 	kf_cos.c kf_sin.c kf_tan.c \
 	ef_acos.c ef_acosh.c ef_asin.c ef_atan2.c \
 	ef_atanh.c ef_cosh.c ef_exp.c ef_fmod.c \
-	erf_gamma.c ef_hypot.c ef_j0.c \
+	ef_tgamma.c ef_hypot.c ef_j0.c \
 	ef_j1.c ef_jn.c erf_lgamma.c \
 	ef_log.c ef_log10.c ef_pow.c ef_rem_pio2.c ef_remainder.c \
 	ef_scalb.c ef_sinh.c ef_sqrt.c \
diff --git a/newlib/libm/math/Makefile.in b/newlib/libm/math/Makefile.in
index 0ac2b1668..7170bf53b 100644
--- a/newlib/libm/math/Makefile.in
+++ b/newlib/libm/math/Makefile.in
@@ -78,7 +78,7 @@ am__objects_1 = lib_a-k_standard.$(OBJEXT) lib_a-k_rem_pio2.$(OBJEXT) \
 	lib_a-e_acosh.$(OBJEXT) lib_a-e_asin.$(OBJEXT) \
 	lib_a-e_atan2.$(OBJEXT) lib_a-e_atanh.$(OBJEXT) \
 	lib_a-e_cosh.$(OBJEXT) lib_a-e_exp.$(OBJEXT) \
-	lib_a-e_fmod.$(OBJEXT) lib_a-er_gamma.$(OBJEXT) \
+	lib_a-e_fmod.$(OBJEXT) lib_a-e_tgamma.$(OBJEXT) \
 	lib_a-e_hypot.$(OBJEXT) lib_a-e_j0.$(OBJEXT) \
 	lib_a-e_j1.$(OBJEXT) lib_a-e_jn.$(OBJEXT) \
 	lib_a-er_lgamma.$(OBJEXT) lib_a-e_log.$(OBJEXT) \
@@ -111,7 +111,7 @@ am__objects_2 = lib_a-kf_rem_pio2.$(OBJEXT) lib_a-kf_cos.$(OBJEXT) \
 	lib_a-ef_asin.$(OBJEXT) lib_a-ef_atan2.$(OBJEXT) \
 	lib_a-ef_atanh.$(OBJEXT) lib_a-ef_cosh.$(OBJEXT) \
 	lib_a-ef_exp.$(OBJEXT) lib_a-ef_fmod.$(OBJEXT) \
-	lib_a-erf_gamma.$(OBJEXT) lib_a-ef_hypot.$(OBJEXT) \
+	lib_a-ef_tgamma.$(OBJEXT) lib_a-ef_hypot.$(OBJEXT) \
 	lib_a-ef_j0.$(OBJEXT) lib_a-ef_j1.$(OBJEXT) \
 	lib_a-ef_jn.$(OBJEXT) lib_a-erf_lgamma.$(OBJEXT) \
 	lib_a-ef_log.$(OBJEXT) lib_a-ef_log10.$(OBJEXT) \
@@ -147,7 +147,7 @@ LTLIBRARIES = $(noinst_LTLIBRARIES)
 libmath_la_LIBADD =
 am__objects_4 = k_standard.lo k_rem_pio2.lo k_cos.lo k_sin.lo k_tan.lo \
 	e_acos.lo e_acosh.lo e_asin.lo e_atan2.lo e_atanh.lo e_cosh.lo \
-	e_exp.lo e_fmod.lo er_gamma.lo e_hypot.lo e_j0.lo e_j1.lo \
+	e_exp.lo e_fmod.lo e_tgamma.lo e_hypot.lo e_j0.lo e_j1.lo \
 	e_jn.lo er_lgamma.lo e_log.lo e_log10.lo e_pow.lo \
 	e_rem_pio2.lo e_remainder.lo e_scalb.lo e_sinh.lo e_sqrt.lo \
 	w_acos.lo w_acosh.lo w_asin.lo w_atan2.lo w_atanh.lo w_cosh.lo \
@@ -159,7 +159,7 @@ am__objects_4 = k_standard.lo k_rem_pio2.lo k_cos.lo k_sin.lo k_tan.lo \
 	s_signif.lo s_sin.lo s_tan.lo s_tanh.lo w_exp2.lo w_tgamma.lo
 am__objects_5 = kf_rem_pio2.lo kf_cos.lo kf_sin.lo kf_tan.lo \
 	ef_acos.lo ef_acosh.lo ef_asin.lo ef_atan2.lo ef_atanh.lo \
-	ef_cosh.lo ef_exp.lo ef_fmod.lo erf_gamma.lo ef_hypot.lo \
+	ef_cosh.lo ef_exp.lo ef_fmod.lo ef_tgamma.lo ef_hypot.lo \
 	ef_j0.lo ef_j1.lo ef_jn.lo erf_lgamma.lo ef_log.lo ef_log10.lo \
 	ef_pow.lo ef_rem_pio2.lo ef_remainder.lo ef_scalb.lo \
 	ef_sinh.lo ef_sqrt.lo wf_acos.lo wf_acosh.lo wf_asin.lo \
@@ -333,7 +333,7 @@ src = k_standard.c k_rem_pio2.c \
 	k_cos.c k_sin.c k_tan.c \
 	e_acos.c e_acosh.c e_asin.c e_atan2.c \
 	e_atanh.c e_cosh.c e_exp.c e_fmod.c \
-	er_gamma.c e_hypot.c e_j0.c \
+	e_tgamma.c e_hypot.c e_j0.c \
 	e_j1.c e_jn.c er_lgamma.c \
 	e_log.c e_log10.c e_pow.c e_rem_pio2.c e_remainder.c \
 	e_scalb.c e_sinh.c e_sqrt.c \
@@ -356,7 +356,7 @@ fsrc = kf_rem_pio2.c \
 	kf_cos.c kf_sin.c kf_tan.c \
 	ef_acos.c ef_acosh.c ef_asin.c ef_atan2.c \
 	ef_atanh.c ef_cosh.c ef_exp.c ef_fmod.c \
-	erf_gamma.c ef_hypot.c ef_j0.c \
+	ef_tgamma.c ef_hypot.c ef_j0.c \
 	ef_j1.c ef_jn.c erf_lgamma.c \
 	ef_log.c ef_log10.c ef_pow.c ef_rem_pio2.c ef_remainder.c \
 	ef_scalb.c ef_sinh.c ef_sqrt.c \
@@ -552,11 +552,11 @@ lib_a-e_fmod.o: e_fmod.c
 lib_a-e_fmod.obj: e_fmod.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_fmod.obj `if test -f 'e_fmod.c'; then $(CYGPATH_W) 'e_fmod.c'; else $(CYGPATH_W) '$(srcdir)/e_fmod.c'; fi`
 
-lib_a-er_gamma.o: er_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-er_gamma.o `test -f 'er_gamma.c' || echo '$(srcdir)/'`er_gamma.c
+lib_a-e_tgamma.o: e_tgamma.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_tgamma.o `test -f 'e_tgamma.c' || echo '$(srcdir)/'`e_tgamma.c
 
-lib_a-er_gamma.obj: er_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-er_gamma.obj `if test -f 'er_gamma.c'; then $(CYGPATH_W) 'er_gamma.c'; else $(CYGPATH_W) '$(srcdir)/er_gamma.c'; fi`
+lib_a-e_tgamma.obj: e_tgamma.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_tgamma.obj `if test -f 'e_tgamma.c'; then $(CYGPATH_W) 'e_tgamma.c'; else $(CYGPATH_W) '$(srcdir)/e_tgamma.c'; fi`
 
 lib_a-e_hypot.o: e_hypot.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-e_hypot.o `test -f 'e_hypot.c' || echo '$(srcdir)/'`e_hypot.c
@@ -948,11 +948,11 @@ lib_a-ef_fmod.o: ef_fmod.c
 lib_a-ef_fmod.obj: ef_fmod.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_fmod.obj `if test -f 'ef_fmod.c'; then $(CYGPATH_W) 'ef_fmod.c'; else $(CYGPATH_W) '$(srcdir)/ef_fmod.c'; fi`
 
-lib_a-erf_gamma.o: erf_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-erf_gamma.o `test -f 'erf_gamma.c' || echo '$(srcdir)/'`erf_gamma.c
+lib_a-ef_tgamma.o: ef_tgamma.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_tgamma.o `test -f 'ef_tgamma.c' || echo '$(srcdir)/'`ef_tgamma.c
 
-lib_a-erf_gamma.obj: erf_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-erf_gamma.obj `if test -f 'erf_gamma.c'; then $(CYGPATH_W) 'erf_gamma.c'; else $(CYGPATH_W) '$(srcdir)/erf_gamma.c'; fi`
+lib_a-ef_tgamma.obj: ef_tgamma.c
+	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_tgamma.obj `if test -f 'ef_tgamma.c'; then $(CYGPATH_W) 'ef_tgamma.c'; else $(CYGPATH_W) '$(srcdir)/ef_tgamma.c'; fi`
 
 lib_a-ef_hypot.o: ef_hypot.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-ef_hypot.o `test -f 'ef_hypot.c' || echo '$(srcdir)/'`ef_hypot.c
diff --git a/newlib/libm/math/er_gamma.c b/newlib/libm/math/e_tgamma.c
index 3c0e241e5..7d8533af8 100644
--- a/newlib/libm/math/er_gamma.c
+++ b/newlib/libm/math/e_tgamma.c
@@ -1,5 +1,5 @@
 
-/* @(#)er_gamma.c 5.1 93/09/24 */
+/* @(#)e_tgamma.c 5.1 93/09/24 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -12,9 +12,8 @@
  *
  */
 
-/* __ieee754_gamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
+/* __ieee754_tgamma(x)
+ * Gamma function. Returns gamma(x)
  *
  * Method: See __ieee754_lgamma_r
  */
@@ -22,11 +21,15 @@
 #include "fdlibm.h"
 
 #ifdef __STDC__
-	double __ieee754_gamma_r(double x, int *signgamp)
+	double __ieee754_tgamma(double x)
 #else
-	double __ieee754_gamma_r(x,signgamp)
-	double x; int *signgamp;
+	double __ieee754_tgamma(x)
+	double x;
 #endif
 {
-	return __ieee754_exp (__ieee754_lgamma_r(x,signgamp));
+	int signgam_local;
+	double y = __ieee754_exp(__ieee754_lgamma_r(x, &signgam_local));
+	if (signgam_local < 0)
+		y = -y;
+	return y;
 }
diff --git a/newlib/libm/math/erf_gamma.c b/newlib/libm/math/ef_tgamma.c
index 9e529dce0..ee429977c 100644
--- a/newlib/libm/math/erf_gamma.c
+++ b/newlib/libm/math/ef_tgamma.c
@@ -1,4 +1,4 @@
-/* erf_gamma.c -- float version of er_gamma.c.
+/* ef_tgamma.c -- float version of e_tgamma.c.
  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
  */
 
@@ -14,9 +14,8 @@
  *
  */
 
-/* __ieee754_gammaf_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
+/* __ieee754_tgammaf(x)
+ * Float version the Gamma function. Returns gamma(x)
  *
  * Method: See __ieee754_lgammaf_r
  */
@@ -24,11 +23,14 @@
 #include "fdlibm.h"
 
 #ifdef __STDC__
-	float __ieee754_gammaf_r(float x, int *signgamp)
+	float __ieee754_tgammaf(float x)
 #else
-	float __ieee754_gammaf_r(x,signgamp)
-	float x; int *signgamp;
+	float __ieee754_tgammaf(x)
 #endif
 {
-	return __ieee754_expf (__ieee754_lgammaf_r(x,signgamp));
+	int signgam_local;
+	float y = __ieee754_expf(__ieee754_lgammaf_r(x, &signgam_local));
+	if (signgam_local < 0)
+		y = -y;
+	return y;
 }
diff --git a/newlib/libm/math/w_gamma.c b/newlib/libm/math/w_gamma.c
index b65d5cc4b..16962c570 100644
--- a/newlib/libm/math/w_gamma.c
+++ b/newlib/libm/math/w_gamma.c
@@ -132,7 +132,7 @@ in terms of the base return values, although the <[signgam]> global for
 /* double gamma(double x)
  * Return the logarithm of the Gamma function of x.
  *
- * Method: call gamma_r
+ * Method: call lgamma
  */
 
 #include "fdlibm.h"
@@ -148,18 +148,7 @@ in terms of the base return values, although the <[signgam]> global for
 	double x;
 #endif
 {
-#ifdef _IEEE_LIBM
-	return __ieee754_gamma_r(x,&(_REENT_SIGNGAM(_REENT)));
-#else
-        double y;
-        y = __ieee754_gamma_r(x,&(_REENT_SIGNGAM(_REENT)));
-        if(_LIB_VERSION == _IEEE_) return y;
-        if(!finite(y)&&finite(x)) {
-	    /* gamma(finite) overflow */
-	    errno = ERANGE;
-        }
-	return y;
-#endif
-}             
+	return lgamma(x);
+}
 
 #endif /* defined(_DOUBLE_IS_32BITS) */
diff --git a/newlib/libm/math/w_tgamma.c b/newlib/libm/math/w_tgamma.c
index c0c011dd0..c09092510 100644
--- a/newlib/libm/math/w_tgamma.c
+++ b/newlib/libm/math/w_tgamma.c
@@ -27,9 +27,7 @@
 #endif
 {
         double y;
-	int local_signgam;
-	y = __ieee754_gamma_r(x,&local_signgam);
-	if (local_signgam < 0) y = -y;
+	y = __ieee754_tgamma(x);
 #ifdef _IEEE_LIBM
 	return y;
 #else
diff --git a/newlib/libm/math/wf_gamma.c b/newlib/libm/math/wf_gamma.c
index f0284a282..29f88a50b 100644
--- a/newlib/libm/math/wf_gamma.c
+++ b/newlib/libm/math/wf_gamma.c
@@ -25,19 +25,8 @@
 	float x;
 #endif
 {
-#ifdef _IEEE_LIBM
-	return __ieee754_gammaf_r(x,&(_REENT_SIGNGAM(_REENT)));
-#else
-        float y;
-        y = __ieee754_gammaf_r(x,&(_REENT_SIGNGAM(_REENT)));
-        if(_LIB_VERSION == _IEEE_) return y;
-        if(!finitef(y)&&finitef(x)) {
-	    /* gammaf(finite) overflow */
-	    errno = ERANGE;
-        }
-	return y;
-#endif
-}             
+	return lgammaf(x);
+}
 
 #ifdef _DOUBLE_IS_32BITS
 
@@ -48,7 +37,7 @@
 	double x;
 #endif
 {
-	return (double) gammaf((float) x);
+	return (double) lgammaf((float) x);
 }
 
 #endif /* defined(_DOUBLE_IS_32BITS) */
diff --git a/newlib/libm/math/wf_tgamma.c b/newlib/libm/math/wf_tgamma.c
index 92df39648..88567b2eb 100644
--- a/newlib/libm/math/wf_tgamma.c
+++ b/newlib/libm/math/wf_tgamma.c
@@ -24,9 +24,7 @@
 #endif
 {
         float y;
-	int local_signgam;
-	y = __ieee754_gammaf_r(x,&local_signgam);
-	if (local_signgam < 0) y = -y;
+	y = __ieee754_tgammaf(x);
 #ifdef _IEEE_LIBM
 	return y;
 #else
diff --git a/newlib/libm/math/wr_gamma.c b/newlib/libm/math/wr_gamma.c
index c4c2a829e..62c2272e9 100644
--- a/newlib/libm/math/wr_gamma.c
+++ b/newlib/libm/math/wr_gamma.c
@@ -27,23 +27,7 @@
 	double x; int *signgamp;
 #endif
 {
-#ifdef _IEEE_LIBM
-	return __ieee754_gamma_r(x,signgamp);
-#else
-        double y;
-        y = __ieee754_gamma_r(x,signgamp);
-        if(_LIB_VERSION == _IEEE_) return y;
-        if(!finite(y)&&finite(x)) {
-	    if(floor(x)==x&&x<=0.0)
-	      /* gamma(-integer) or gamma(0) */
-	      errno = EDOM;
-	    else
-	      /* gamma(finite) overflow */
-	      errno = ERANGE;
-	    return HUGE_VALF;
-        } else
-            return y;
-#endif
+	return lgamma_r(x, signgamp);
 }
 
 #endif /* defined(_DOUBLE_IS_32BITS) */
diff --git a/newlib/libm/math/wr_lgamma.c b/newlib/libm/math/wr_lgamma.c
index d961db58d..7f86f37fb 100644
--- a/newlib/libm/math/wr_lgamma.c
+++ b/newlib/libm/math/wr_lgamma.c
@@ -40,9 +40,8 @@
 	    else
 	      /* lgamma(finite) overflow */
 	      errno = ERANGE;
-	    return HUGE_VAL;
-        } else
-            return y;
+        }
+	return y;
 #endif
 }
 
diff --git a/newlib/libm/math/wrf_gamma.c b/newlib/libm/math/wrf_gamma.c
index d43c7f03d..15ec1eee9 100644
--- a/newlib/libm/math/wrf_gamma.c
+++ b/newlib/libm/math/wrf_gamma.c
@@ -27,22 +27,5 @@
 	float x; int *signgamp;
 #endif
 {
-#ifdef _IEEE_LIBM
-	return __ieee754_gammaf_r(x,signgamp);
-#else
-        float y;
-        y = __ieee754_gammaf_r(x,signgamp);
-        if(_LIB_VERSION == _IEEE_) return y;
-        if(!finitef(y)&&finitef(x)) {
-	    if(floorf(x)==x&&x<=0.0f) {
-		/* gammaf(-integer) or gamma(0) */
-		errno = EDOM;
-	    } else {
-		/* gammaf(finite) overflow */
-		errno = ERANGE;
-	    }
-	    return HUGE_VALF;
-        } else
-            return y;
-#endif
-}             
+	return lgammaf_r(x, signgamp);
+}
diff --git a/newlib/libm/math/wrf_lgamma.c b/newlib/libm/math/wrf_lgamma.c
index 2785b87e5..8822cdf78 100644
--- a/newlib/libm/math/wrf_lgamma.c
+++ b/newlib/libm/math/wrf_lgamma.c
@@ -41,8 +41,7 @@
 		/* lgammaf(finite) overflow */
 		errno = ERANGE;
 	    }
-	    return HUGE_VALF;
-        } else
-            return y;
+        }
+	return y;
 #endif
-}             
+}