[PATCH 0/3] libm: Clean up gamma functions

[PATCH 0/3] libm: Clean up gamma functions

[Keith Packard]

This  series  cleans   up  some  API  oddities   and  inaccuracies  in
errno/exception reporting  from the  various gamma  related functions.
With this series applied, newlib now  matches the POSIX spec (and glibc)
for  tgamma/lgamma  functions. I've  left  the  'gamma' functions  BSD
compatible (tgamma instead of lgamma) to avoid changing the API.

[PATCH 1/3] libm: Fix sign value returned from __ieee754_lgamma*_r(-0)

	Following the C spec, gamma(-0) should be -INFINITY rather than
	+INFINITY, so the *signgamp value in lgamma should be -1
	rather than +1 for this case.

[PATCH 2/3] libm: Remove __ieee754_gamma_r variants

	Back in 2002, a patch changed the various gamma functions so
	that the returned abs(gamma) instead of the ln of gamma. This
	left the _r variants in place, with their extra parameter to
	return the sign of the resulting value. This doesn't make any
	sense to me and looks like it was just a mistake. I've changed
	the __ieee754_gamma functions to apply the sign value locally
	and eliminated the gamma_r versions as those are no longer useful.

[PATCH 3/3] libm: Adjust errno/exception values for gamma/lgamma

	There's a pretty detailed discussion in the commit message
	here, the short version is that tgamma and lgamma are
	specified with different errno/exception behaviour for
	one class of inputs (negative integers).

[PATCH 1/3] libm: Fix sign value returned from __ieee754_lgamma*_r(-0)

[Keith Packard]

The sign of the INFINITY returned from these cases needs to match the
sign of the zero.

Signed-off-by: Keith Packard <keithp@keithp.com>
---
 newlib/libm/math/er_lgamma.c  | 6 +++++-
 newlib/libm/math/erf_lgamma.c | 6 +++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/newlib/libm/math/er_lgamma.c b/newlib/libm/math/er_lgamma.c
index 386a8a73b..36408382f 100644
--- a/newlib/libm/math/er_lgamma.c
+++ b/newlib/libm/math/er_lgamma.c
@@ -225,7 +225,11 @@ static double zero=  0.00000000000000000000e+00;
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7ff00000) return x*x;
-	if((ix|lx)==0) return one/zero;
+	if((ix|lx)==0) {
+	    if(hx<0)
+	        *signgamp = -1;
+	    return one/zero;
+	}
 	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
index 3c6ba02af..a45423949 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -160,7 +160,11 @@ static float zero=  0.0000000000e+00;
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7f800000) return x*x;
-	if(ix==0) return one/zero;
+	if(ix==0) {
+	    if(hx<0)
+	        *signgamp = -1;
+	    return one/zero;
+	}
 	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
-- 
2.28.0

[PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

Gamma should consume the sign reported by lgamma internally instead of
providing it back to the application. This removes the need to have _r
variants as there is no longer any re-entrancy issues with the API.

Signed-off-by: Keith Packard <keithp@keithp.com>
---
 newlib/libc/include/math.h                 |  2 -
 newlib/libc/sys/linux/cmath/math_private.h |  2 -
 newlib/libm/common/fdlibm.h                |  4 +-
 newlib/libm/math/Makefile.am               |  4 +-
 newlib/libm/math/Makefile.in               | 24 +++--------
 newlib/libm/math/er_gamma.c                | 12 ++++--
 newlib/libm/math/erf_gamma.c               | 12 ++++--
 newlib/libm/math/w_gamma.c                 | 15 +------
 newlib/libm/math/w_tgamma.c                |  4 +-
 newlib/libm/math/wf_gamma.c                | 15 +------
 newlib/libm/math/wf_tgamma.c               |  4 +-
 newlib/libm/math/wr_gamma.c                | 49 ----------------------
 newlib/libm/math/wrf_gamma.c               | 48 ---------------------
 13 files changed, 32 insertions(+), 163 deletions(-)
 delete mode 100644 newlib/libm/math/wr_gamma.c
 delete mode 100644 newlib/libm/math/wrf_gamma.c

diff --git a/newlib/libc/include/math.h b/newlib/libc/include/math.h
index 3399d3649..28d4df10f 100644
--- a/newlib/libc/include/math.h
+++ b/newlib/libc/include/math.h
@@ -515,9 +515,7 @@ extern float dremf (float, float);
 #ifdef __CYGWIN__
 extern float dreml (long double, long double);
 #endif /* __CYGWIN__ */
-extern double gamma_r (double, int *);
 extern double lgamma_r (double, int *);
-extern float gammaf_r (float, int *);
 extern float lgammaf_r (float, int *);
 #endif
 
diff --git a/newlib/libc/sys/linux/cmath/math_private.h b/newlib/libc/sys/linux/cmath/math_private.h
index 3e32b29ba..f5792f7fa 100644
--- a/newlib/libc/sys/linux/cmath/math_private.h
+++ b/newlib/libc/sys/linux/cmath/math_private.h
@@ -178,7 +178,6 @@ extern double __ieee754_cosh (double);
 extern double __ieee754_fmod (double,double);
 extern double __ieee754_pow (double,double);
 extern double __ieee754_lgamma_r (double,int *);
-extern double __ieee754_gamma_r (double,int *);
 extern double __ieee754_lgamma (double);
 extern double __ieee754_gamma (double);
 extern double __ieee754_log10 (double);
@@ -241,7 +240,6 @@ extern float __ieee754_coshf (float);
 extern float __ieee754_fmodf (float,float);
 extern float __ieee754_powf (float,float);
 extern float __ieee754_lgammaf_r (float,int *);
-extern float __ieee754_gammaf_r (float,int *);
 extern float __ieee754_lgammaf (float);
 extern float __ieee754_gammaf (float);
 extern float __ieee754_log10f (float);
diff --git a/newlib/libm/common/fdlibm.h b/newlib/libm/common/fdlibm.h
index 8dffc832d..5226c2e13 100644
--- a/newlib/libm/common/fdlibm.h
+++ b/newlib/libm/common/fdlibm.h
@@ -159,7 +159,7 @@ extern double __ieee754_cosh __P((double));
 extern double __ieee754_fmod __P((double,double));
 extern double __ieee754_pow __P((double,double));
 extern double __ieee754_lgamma_r __P((double,int *));
-extern double __ieee754_gamma_r __P((double,int *));
+extern double __ieee754_gamma __P((double));
 extern double __ieee754_log10 __P((double));
 extern double __ieee754_sinh __P((double));
 extern double __ieee754_hypot __P((double,double));
@@ -205,7 +205,7 @@ extern float __ieee754_coshf __P((float));
 extern float __ieee754_fmodf __P((float,float));
 extern float __ieee754_powf __P((float,float));
 extern float __ieee754_lgammaf_r __P((float,int *));
-extern float __ieee754_gammaf_r __P((float,int *));
+extern float __ieee754_gammaf __P((float));
 extern float __ieee754_log10f __P((float));
 extern float __ieee754_sinhf __P((float));
 extern float __ieee754_hypotf __P((float,float));
diff --git a/newlib/libm/math/Makefile.am b/newlib/libm/math/Makefile.am
index e745159ae..87e47dbd5 100644
--- a/newlib/libm/math/Makefile.am
+++ b/newlib/libm/math/Makefile.am
@@ -14,7 +14,7 @@ src = 	k_standard.c k_rem_pio2.c \
 	e_scalb.c e_sinh.c e_sqrt.c \
 	w_acos.c w_acosh.c w_asin.c w_atan2.c \
 	w_atanh.c w_cosh.c w_exp.c w_fmod.c \
-	w_gamma.c wr_gamma.c w_hypot.c w_j0.c \
+	w_gamma.c w_hypot.c w_j0.c \
 	w_j1.c w_jn.c w_lgamma.c wr_lgamma.c \
 	w_log.c w_log10.c w_pow.c w_remainder.c \
 	w_scalb.c w_sinh.c w_sqrt.c \
@@ -37,7 +37,7 @@ fsrc =	kf_rem_pio2.c \
 	ef_scalb.c ef_sinh.c ef_sqrt.c \
 	wf_acos.c wf_acosh.c wf_asin.c wf_atan2.c \
 	wf_atanh.c wf_cosh.c wf_exp.c wf_fmod.c \
-	wf_gamma.c wrf_gamma.c wf_hypot.c wf_j0.c \
+	wf_gamma.c wf_hypot.c wf_j0.c \
 	wf_j1.c wf_jn.c wf_lgamma.c wrf_lgamma.c \
 	wf_log.c wf_log10.c wf_pow.c wf_remainder.c \
 	wf_scalb.c wf_sinh.c wf_sqrt.c \
diff --git a/newlib/libm/math/Makefile.in b/newlib/libm/math/Makefile.in
index 0ac2b1668..9ba1be8ca 100644
--- a/newlib/libm/math/Makefile.in
+++ b/newlib/libm/math/Makefile.in
@@ -90,7 +90,7 @@ am__objects_1 = lib_a-k_standard.$(OBJEXT) lib_a-k_rem_pio2.$(OBJEXT) \
 	lib_a-w_atan2.$(OBJEXT) lib_a-w_atanh.$(OBJEXT) \
 	lib_a-w_cosh.$(OBJEXT) lib_a-w_exp.$(OBJEXT) \
 	lib_a-w_fmod.$(OBJEXT) lib_a-w_gamma.$(OBJEXT) \
-	lib_a-wr_gamma.$(OBJEXT) lib_a-w_hypot.$(OBJEXT) \
+	lib_a-w_hypot.$(OBJEXT) \
 	lib_a-w_j0.$(OBJEXT) lib_a-w_j1.$(OBJEXT) lib_a-w_jn.$(OBJEXT) \
 	lib_a-w_lgamma.$(OBJEXT) lib_a-wr_lgamma.$(OBJEXT) \
 	lib_a-w_log.$(OBJEXT) lib_a-w_log10.$(OBJEXT) \
@@ -122,7 +122,7 @@ am__objects_2 = lib_a-kf_rem_pio2.$(OBJEXT) lib_a-kf_cos.$(OBJEXT) \
 	lib_a-wf_asin.$(OBJEXT) lib_a-wf_atan2.$(OBJEXT) \
 	lib_a-wf_atanh.$(OBJEXT) lib_a-wf_cosh.$(OBJEXT) \
 	lib_a-wf_exp.$(OBJEXT) lib_a-wf_fmod.$(OBJEXT) \
-	lib_a-wf_gamma.$(OBJEXT) lib_a-wrf_gamma.$(OBJEXT) \
+	lib_a-wf_gamma.$(OBJEXT) \
 	lib_a-wf_hypot.$(OBJEXT) lib_a-wf_j0.$(OBJEXT) \
 	lib_a-wf_j1.$(OBJEXT) lib_a-wf_jn.$(OBJEXT) \
 	lib_a-wf_lgamma.$(OBJEXT) lib_a-wrf_lgamma.$(OBJEXT) \
@@ -151,7 +151,7 @@ am__objects_4 = k_standard.lo k_rem_pio2.lo k_cos.lo k_sin.lo k_tan.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 \
-	w_exp.lo w_fmod.lo w_gamma.lo wr_gamma.lo w_hypot.lo w_j0.lo \
+	w_exp.lo w_fmod.lo w_gamma.lo w_hypot.lo w_j0.lo \
 	w_j1.lo w_jn.lo w_lgamma.lo wr_lgamma.lo w_log.lo w_log10.lo \
 	w_pow.lo w_remainder.lo w_scalb.lo w_sinh.lo w_sqrt.lo \
 	w_sincos.lo w_drem.lo s_asinh.lo s_atan.lo s_ceil.lo s_cos.lo \
@@ -164,7 +164,7 @@ am__objects_5 = kf_rem_pio2.lo kf_cos.lo kf_sin.lo kf_tan.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 \
 	wf_atan2.lo wf_atanh.lo wf_cosh.lo wf_exp.lo wf_fmod.lo \
-	wf_gamma.lo wrf_gamma.lo wf_hypot.lo wf_j0.lo wf_j1.lo \
+	wf_gamma.lo wf_hypot.lo wf_j0.lo wf_j1.lo \
 	wf_jn.lo wf_lgamma.lo wrf_lgamma.lo wf_log.lo wf_log10.lo \
 	wf_pow.lo wf_remainder.lo wf_scalb.lo wf_sinh.lo wf_sqrt.lo \
 	wf_sincos.lo wf_drem.lo sf_asinh.lo sf_atan.lo sf_ceil.lo \
@@ -339,7 +339,7 @@ src = k_standard.c k_rem_pio2.c \
 	e_scalb.c e_sinh.c e_sqrt.c \
 	w_acos.c w_acosh.c w_asin.c w_atan2.c \
 	w_atanh.c w_cosh.c w_exp.c w_fmod.c \
-	w_gamma.c wr_gamma.c w_hypot.c w_j0.c \
+	w_gamma.c w_hypot.c w_j0.c \
 	w_j1.c w_jn.c w_lgamma.c wr_lgamma.c \
 	w_log.c w_log10.c w_pow.c w_remainder.c \
 	w_scalb.c w_sinh.c w_sqrt.c \
@@ -362,7 +362,7 @@ fsrc = kf_rem_pio2.c \
 	ef_scalb.c ef_sinh.c ef_sqrt.c \
 	wf_acos.c wf_acosh.c wf_asin.c wf_atan2.c \
 	wf_atanh.c wf_cosh.c wf_exp.c wf_fmod.c \
-	wf_gamma.c wrf_gamma.c wf_hypot.c wf_j0.c \
+	wf_gamma.c wf_hypot.c wf_j0.c \
 	wf_j1.c wf_jn.c wf_lgamma.c wrf_lgamma.c \
 	wf_log.c wf_log10.c wf_pow.c wf_remainder.c \
 	wf_scalb.c wf_sinh.c wf_sqrt.c \
@@ -690,12 +690,6 @@ lib_a-w_gamma.o: w_gamma.c
 lib_a-w_gamma.obj: w_gamma.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-w_gamma.obj `if test -f 'w_gamma.c'; then $(CYGPATH_W) 'w_gamma.c'; else $(CYGPATH_W) '$(srcdir)/w_gamma.c'; fi`
 
-lib_a-wr_gamma.o: wr_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-wr_gamma.o `test -f 'wr_gamma.c' || echo '$(srcdir)/'`wr_gamma.c
-
-lib_a-wr_gamma.obj: wr_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-wr_gamma.obj `if test -f 'wr_gamma.c'; then $(CYGPATH_W) 'wr_gamma.c'; else $(CYGPATH_W) '$(srcdir)/wr_gamma.c'; fi`
-
 lib_a-w_hypot.o: w_hypot.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-w_hypot.o `test -f 'w_hypot.c' || echo '$(srcdir)/'`w_hypot.c
 
@@ -1086,12 +1080,6 @@ lib_a-wf_gamma.o: wf_gamma.c
 lib_a-wf_gamma.obj: wf_gamma.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-wf_gamma.obj `if test -f 'wf_gamma.c'; then $(CYGPATH_W) 'wf_gamma.c'; else $(CYGPATH_W) '$(srcdir)/wf_gamma.c'; fi`
 
-lib_a-wrf_gamma.o: wrf_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-wrf_gamma.o `test -f 'wrf_gamma.c' || echo '$(srcdir)/'`wrf_gamma.c
-
-lib_a-wrf_gamma.obj: wrf_gamma.c
-	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-wrf_gamma.obj `if test -f 'wrf_gamma.c'; then $(CYGPATH_W) 'wrf_gamma.c'; else $(CYGPATH_W) '$(srcdir)/wrf_gamma.c'; fi`
-
 lib_a-wf_hypot.o: wf_hypot.c
 	$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-wf_hypot.o `test -f 'wf_hypot.c' || echo '$(srcdir)/'`wf_hypot.c
 
diff --git a/newlib/libm/math/er_gamma.c b/newlib/libm/math/er_gamma.c
index 3c0e241e5..9252e2d04 100644
--- a/newlib/libm/math/er_gamma.c
+++ b/newlib/libm/math/er_gamma.c
@@ -22,11 +22,15 @@
 #include "fdlibm.h"
 
 #ifdef __STDC__
-	double __ieee754_gamma_r(double x, int *signgamp)
+	double __ieee754_gamma(double x)
 #else
-	double __ieee754_gamma_r(x,signgamp)
-	double x; int *signgamp;
+	double __ieee754_gamma(x)
+	double x;
 #endif
 {
-	return __ieee754_exp (__ieee754_lgamma_r(x,signgamp));
+	int local_signgam = 1;
+	double y = __ieee754_exp (__ieee754_lgamma_r(x, &local_signgam));
+	if (local_signgam < 0)
+		y = -y;
+	return y;
 }
diff --git a/newlib/libm/math/erf_gamma.c b/newlib/libm/math/erf_gamma.c
index 9e529dce0..bddac60c7 100644
--- a/newlib/libm/math/erf_gamma.c
+++ b/newlib/libm/math/erf_gamma.c
@@ -24,11 +24,15 @@
 #include "fdlibm.h"
 
 #ifdef __STDC__
-	float __ieee754_gammaf_r(float x, int *signgamp)
+	float __ieee754_gammaf(float x)
 #else
-	float __ieee754_gammaf_r(x,signgamp)
-	float x; int *signgamp;
+	float __ieee754_gammaf(x)
+	float x;
 #endif
 {
-	return __ieee754_expf (__ieee754_lgammaf_r(x,signgamp));
+	int local_signgam = 1;
+	float y = __ieee754_expf (__ieee754_lgammaf_r(x,&local_signgam));
+	if (local_signgam < 0)
+		y = -y;
+	return y;
 }
diff --git a/newlib/libm/math/w_gamma.c b/newlib/libm/math/w_gamma.c
index b65d5cc4b..bb47fab31 100644
--- a/newlib/libm/math/w_gamma.c
+++ b/newlib/libm/math/w_gamma.c
@@ -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 tgamma(x);
+}
 
 #endif /* defined(_DOUBLE_IS_32BITS) */
diff --git a/newlib/libm/math/w_tgamma.c b/newlib/libm/math/w_tgamma.c
index c0c011dd0..cbb3e7f7f 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_gamma(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..4e961c864 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 tgammaf(x);
+}
 
 #ifdef _DOUBLE_IS_32BITS
 
diff --git a/newlib/libm/math/wf_tgamma.c b/newlib/libm/math/wf_tgamma.c
index 92df39648..8f0c501b8 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_gammaf(x);
 #ifdef _IEEE_LIBM
 	return y;
 #else
diff --git a/newlib/libm/math/wr_gamma.c b/newlib/libm/math/wr_gamma.c
deleted file mode 100644
index c4c2a829e..000000000
--- a/newlib/libm/math/wr_gamma.c
+++ /dev/null
@@ -1,49 +0,0 @@
-
-/* @(#)wr_gamma.c 5.1 93/09/24 */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
- * is preserved.
- * ====================================================
- */
-
-/* 
- * wrapper double gamma_r(double x, int *signgamp)
- */
-
-#include "fdlibm.h"
-#include <errno.h>
-
-#ifndef _DOUBLE_IS_32BITS
-
-#ifdef __STDC__
-	double gamma_r(double x, int *signgamp) /* wrapper lgamma_r */
-#else
-	double gamma_r(x,signgamp)              /* wrapper lgamma_r */
-	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
-}
-
-#endif /* defined(_DOUBLE_IS_32BITS) */
diff --git a/newlib/libm/math/wrf_gamma.c b/newlib/libm/math/wrf_gamma.c
deleted file mode 100644
index d43c7f03d..000000000
--- a/newlib/libm/math/wrf_gamma.c
+++ /dev/null
@@ -1,48 +0,0 @@
-/* wrf_gamma.c -- float version of wr_gamma.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
- * is preserved.
- * ====================================================
- */
-
-/* 
- * wrapper float gammaf_r(float x, int *signgamp)
- */
-
-#include "fdlibm.h"
-#include <errno.h>
-
-#ifdef __STDC__
-	float gammaf_r(float x, int *signgamp) /* wrapper lgammaf_r */
-#else
-	float gammaf_r(x,signgamp)              /* wrapper lgammaf_r */
-	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
-}             
-- 
2.28.0

[PATCH 3/3] libm: Adjust errno/exception values for gamma/lgamma

[Keith Packard]

This makes all of gamma/lgamma functions match POSIX (and glibc) for
errno and exception values for both gamma and lgamma functions.

The big change is to support the exception/errno value differences for
tgamma/lgamma for negative integer arguments. tgamma produces
EDOM/FE_INVALID while lgamma produces ERANGE/FE_DIVBYZERO.

This was done by splitting the lowest level __ieee754_lgamma*_r
functions apart. The new lower-level ___ieee754_lgamma*_r functions
can now produce exceptions which are correct for either lgamma or
tgamma and use the value in the signgamp parameter to select which
mode to operate in.

All functions now return EDOM/FE_INVALID for -INFINITY, although POSIX
doesn't specify this behavior for lgamma. It does match glibc at
least.

Signed-off-by: Keith Packard <keithp@keithp.com>
---
 newlib/libm/common/fdlibm.h   |  2 ++
 newlib/libm/math/er_gamma.c   |  2 +-
 newlib/libm/math/er_lgamma.c  | 32 ++++++++++++++++++++++++++------
 newlib/libm/math/erf_gamma.c  |  2 +-
 newlib/libm/math/erf_lgamma.c | 30 +++++++++++++++++++++++++-----
 newlib/libm/math/k_standard.c | 13 ++++++++++---
 newlib/libm/math/w_tgamma.c   |  6 ++++--
 newlib/libm/math/wf_tgamma.c  |  7 +++++--
 8 files changed, 74 insertions(+), 20 deletions(-)

diff --git a/newlib/libm/common/fdlibm.h b/newlib/libm/common/fdlibm.h
index 5226c2e13..5613747bb 100644
--- a/newlib/libm/common/fdlibm.h
+++ b/newlib/libm/common/fdlibm.h
@@ -159,6 +159,7 @@ extern double __ieee754_cosh __P((double));
 extern double __ieee754_fmod __P((double,double));
 extern double __ieee754_pow __P((double,double));
 extern double __ieee754_lgamma_r __P((double,int *));
+extern double ___ieee754_lgamma_r __P((double,int *));
 extern double __ieee754_gamma __P((double));
 extern double __ieee754_log10 __P((double));
 extern double __ieee754_sinh __P((double));
@@ -205,6 +206,7 @@ extern float __ieee754_coshf __P((float));
 extern float __ieee754_fmodf __P((float,float));
 extern float __ieee754_powf __P((float,float));
 extern float __ieee754_lgammaf_r __P((float,int *));
+extern float ___ieee754_lgammaf_r __P((float,int *));
 extern float __ieee754_gammaf __P((float));
 extern float __ieee754_log10f __P((float));
 extern float __ieee754_sinhf __P((float));
diff --git a/newlib/libm/math/er_gamma.c b/newlib/libm/math/er_gamma.c
index 9252e2d04..d7731c1c3 100644
--- a/newlib/libm/math/er_gamma.c
+++ b/newlib/libm/math/er_gamma.c
@@ -29,7 +29,7 @@
 #endif
 {
 	int local_signgam = 1;
-	double y = __ieee754_exp (__ieee754_lgamma_r(x, &local_signgam));
+	double y = __ieee754_exp (___ieee754_lgamma_r(x, &local_signgam));
 	if (local_signgam < 0)
 		y = -y;
 	return y;
diff --git a/newlib/libm/math/er_lgamma.c b/newlib/libm/math/er_lgamma.c
index 36408382f..a93d54cf2 100644
--- a/newlib/libm/math/er_lgamma.c
+++ b/newlib/libm/math/er_lgamma.c
@@ -210,21 +210,26 @@ static double zero=  0.00000000000000000000e+00;
 
 
 #ifdef __STDC__
-	double __ieee754_lgamma_r(double x, int *signgamp)
+	double ___ieee754_lgamma_r(double x, int *signgamp)
 #else
-	double __ieee754_lgamma_r(x,signgamp)
+	double ___ieee754_lgamma_r(x, signgamp)
 	double x; int *signgamp;
 #endif
 {
 	double t,y,z,nadj = 0.0,p,p1,p2,p3,q,r,w;
 	__int32_t i,hx,lx,ix;
+	int mode = *signgamp;
 
 	EXTRACT_WORDS(hx,lx,x);
 
     /* purge off +-inf, NaN, +-0, and negative arguments */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
-	if(ix>=0x7ff00000) return x*x;
+	if(ix>=0x7ff00000) {
+	    if (hx<0 && mode)
+		return __math_invalid(x);
+	    return x*x;
+	}
 	if((ix|lx)==0) {
 	    if(hx<0)
 	        *signgamp = -1;
@@ -237,10 +242,19 @@ static double zero=  0.00000000000000000000e+00;
 	    } else return -__ieee754_log(x);
 	}
 	if(hx<0) {
-	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */
-		return one/zero;
+	    if(ix>=0x43300000) { /* |x|>=2**52, must be -integer */
+		if (mode)
+		    return __math_invalid(x);
+		else
+		    return one/zero; /* -integer */
+	    }
 	    t = sin_pi(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) {
+		if (mode)
+		    return __math_invalid(x);
+		else
+		    return one/zero; /* -integer */
+	    }
 	    nadj = __ieee754_log(pi/fabs(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
@@ -311,3 +325,9 @@ static double zero=  0.00000000000000000000e+00;
 	if(hx<0) r = nadj - r;
 	return r;
 }
+
+double __ieee754_lgamma_r(double x, int *signgamp)
+{
+    *signgamp = 0;
+    return ___ieee754_lgamma_r(x, signgamp);
+}
diff --git a/newlib/libm/math/erf_gamma.c b/newlib/libm/math/erf_gamma.c
index bddac60c7..9eeb5d978 100644
--- a/newlib/libm/math/erf_gamma.c
+++ b/newlib/libm/math/erf_gamma.c
@@ -31,7 +31,7 @@
 #endif
 {
 	int local_signgam = 1;
-	float y = __ieee754_expf (__ieee754_lgammaf_r(x,&local_signgam));
+	float y = __ieee754_expf (___ieee754_lgammaf_r(x,&local_signgam));
 	if (local_signgam < 0)
 		y = -y;
 	return y;
diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
index a45423949..ec83c82ff 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -145,21 +145,26 @@ static float zero=  0.0000000000e+00;
 
 
 #ifdef __STDC__
-	float __ieee754_lgammaf_r(float x, int *signgamp)
+	float ___ieee754_lgammaf_r(float x, int *signgamp)
 #else
-	float __ieee754_lgammaf_r(x,signgamp)
+	float ___ieee754_lgammaf_r(x,signgamp)
 	float x; int *signgamp;
 #endif
 {
 	float t,y,z,nadj = 0.0,p,p1,p2,p3,q,r,w;
 	__int32_t i,hx,ix;
+	int mode = *signgamp;
 
 	GET_FLOAT_WORD(hx,x);
 
     /* purge off +-inf, NaN, +-0, and negative arguments */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) return x*x;
+	if(ix>=0x7f800000) {
+	    if (hx<0 && mode)
+		return __math_invalidf(x);
+	    return x*x;
+	}
 	if(ix==0) {
 	    if(hx<0)
 	        *signgamp = -1;
@@ -172,10 +177,18 @@ static float zero=  0.0000000000e+00;
 	    } else return -__ieee754_logf(x);
 	}
 	if(hx<0) {
-	    if(ix>=0x4b000000) 	/* |x|>=2**23, must be -integer */
+	    if(ix>=0x4b000000) { 	/* |x|>=2**23, must be -integer */
+		if (mode)
+		    return __math_invalidf(x);
 		return one/zero;
+	    }
 	    t = sin_pif(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) {
+		/* tgamma wants NaN instead of INFINITY */
+		if (mode)
+		    return __math_invalidf(x);
+		return one/zero; /* -integer */
+	    }
 	    nadj = __ieee754_logf(pi/fabsf(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
@@ -246,3 +259,10 @@ static float zero=  0.0000000000e+00;
 	if(hx<0) r = nadj - r;
 	return r;
 }
+
+
+float __ieee754_lgammaf_r(float x, int *signgamp)
+{
+    *signgamp = 0;
+    return ___ieee754_lgammaf_r(x, signgamp);
+}
diff --git a/newlib/libm/math/k_standard.c b/newlib/libm/math/k_standard.c
index 906412ba7..8380ad682 100644
--- a/newlib/libm/math/k_standard.c
+++ b/newlib/libm/math/k_standard.c
@@ -74,7 +74,8 @@ static double zero = 0.0;	/* used as const */
  *	39-- yn(x>X_TLOSS, n)
  *	40-- gamma(finite) overflow
  *	41-- gamma(-integer)
- *	42-- pow(NaN,0.0)
+ *	42-- gamma(0)
+ *	43-- pow(NaN,0.0)
  */
 
 
@@ -336,12 +337,18 @@ static double zero = 0.0;	/* used as const */
 		break;
 	    case 41:
 	    case 141:
-		/* gamma(-integer) or gamma(0) */
-		retval = HUGE_VAL;
+		/* gamma(-integer) */
+		retval = (double) NAN;
 		errno = EDOM;
 		break;
 	    case 42:
 	    case 142:
+		/* gamma(0) */
+		retval = copysign(HUGE_VAL,x);
+		errno = ERANGE;
+		break;
+	    case 43:
+	    case 143:
 		/* pow(NaN,0.0) */
 		/* Not an error.  */
 		retval = 1.0;
diff --git a/newlib/libm/math/w_tgamma.c b/newlib/libm/math/w_tgamma.c
index cbb3e7f7f..af8d1e918 100644
--- a/newlib/libm/math/w_tgamma.c
+++ b/newlib/libm/math/w_tgamma.c
@@ -34,8 +34,10 @@
 	if(_LIB_VERSION == _IEEE_) return y;
 
 	if(!finite(y)&&finite(x)) {
-	  if(floor(x)==x&&x<=0.0)
-	    return __kernel_standard(x,x,41); /* tgamma pole */
+	  if (x==0.0)
+	    return __kernel_standard(x,x,42); /* tgamma pole */
+	  else if(floor(x)==x&&x<0.0)
+	    return __kernel_standard(x,x,41); /* tgamma domain */
 	  else
 	    return __kernel_standard(x,x,40); /* tgamma overflow */
 	}
diff --git a/newlib/libm/math/wf_tgamma.c b/newlib/libm/math/wf_tgamma.c
index 8f0c501b8..bf6273b60 100644
--- a/newlib/libm/math/wf_tgamma.c
+++ b/newlib/libm/math/wf_tgamma.c
@@ -31,8 +31,11 @@
 	if(_LIB_VERSION == _IEEE_) return y;
 
 	if(!finitef(y)&&finitef(x)) {
-	  if(floorf(x)==x&&x<=(float)0.0)
-	    /* tgammaf pole */
+	  if (x==0.0f)
+	    /* tgammf pole */
+	    return (float)__kernel_standard((double)x,(double)x,142);
+	  else if(floorf(x)==x&&x<(float)0.0)
+	    /* tgammaf domain */
 	    return (float)__kernel_standard((double)x,(double)x,141);
 	  else
 	    /* tgammaf overflow */
-- 
2.28.0

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Corinna Vinschen]

Hi Keith,

On Aug 26 10:03, Keith Packard via Newlib wrote:
> Gamma should consume the sign reported by lgamma internally instead of
> providing it back to the application. This removes the need to have _r
> variants as there is no longer any re-entrancy issues with the API.

You can't do that.  gamma_r/gammaf_r/lgamma_r/lgammaf_r are BSD
functions.  They have been exported by Cygwin since 2001.  The entry
points need to be kept available with unchanged semantics.


Corinna

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 413 bytes --]

Corinna Vinschen via Newlib <newlib@sourceware.org> writes:

> You can't do that.  gamma_r/gammaf_r/lgamma_r/lgammaf_r are BSD
> functions.  They have been exported by Cygwin since 2001.  The entry
> points need to be kept available with unchanged semantics.

They changed on accident in 2002 -- gamma_r *used* to be an alias for
lgamma_r and was changed to be a (broken) alias for tgamma.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Corinna Vinschen]

On Aug 26 12:10, Keith Packard via Newlib wrote:
> Corinna Vinschen via Newlib <newlib@sourceware.org> writes:
> 
> > You can't do that.  gamma_r/gammaf_r/lgamma_r/lgammaf_r are BSD
> > functions.  They have been exported by Cygwin since 2001.  The entry
> > points need to be kept available with unchanged semantics.
> 
> They changed on accident in 2002 -- gamma_r *used* to be an alias for
> lgamma_r and was changed to be a (broken) alias for tgamma.

Nevertheless, the symbols have to be exported for backward compatibility.
So you're saying aliasing gamma_r to lgamma_r and gammaf_r to lgammaf_r
in the DLLs export table would be sufficient?


Thanks,
Corinna

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 3459 bytes --]

Corinna Vinschen via Newlib <newlib@sourceware.org> writes:

> Nevertheless, the symbols have to be exported for backward compatibility.
> So you're saying aliasing gamma_r to lgamma_r and gammaf_r to lgammaf_r
> in the DLLs export table would be sufficient?

It depends if you want the pre-2002 functionality or the post-2002
functionality.

Before 2002, gamma_r and gammaf_r were as BSD originally specified them,
returning the log of gamma (ln(|Γ(x)|)) with the sign of gamma stored
through the second, int *, argument. In 2002, this patch was made to
er_gamma.c:

$ git show 0953fe640f177b565578ed7ecc77169ec1a914fa er_gamma.c

diff --git a/newlib/libm/math/er_gamma.c b/newlib/libm/math/er_gamma.c
index a7183c50f..3c0e241e5 100644
@@ -28,5 +28,5 @@
 	double x; int *signgamp;
 #endif
 {
-	return __ieee754_lgamma_r(x,signgamp);
+	return __ieee754_exp (__ieee754_lgamma_r(x,signgamp));
 }

Before this patch, __ieee754_gamma_r was simply another name for
__ieee754_lgamma_r. After this patch, __ieee754_gamma_r became something
that doesn't exist in any other math library: it returns |Γ(x)| and
stores the sign through the second, int *, argument. Internally, this is
used in w_gamma.c, w_tgamma.c, wr_gamma.c.

w_tgamma.c uses it correctly:

	y = __ieee754_gamma_r(x,&local_signgam);
	if (local_signgam < 0) y = -y;

w_gamma.c, which exports the 'gamma' function uses it *incorrectly* --
it didn't change after the above patch, and so it unexpectedly changed
from returning ln(|Γ(x)|) to returning |Γ(x)|. Applications using
'gamma' instead of 'tgamma' will be getting the wrong answer for some
parameters.

I have to assume this was just an oversight. In 2009, another patch to
w_gamma.c adds a comment explicitly stating that gamma and gamma_r
return ln(|Γ(x)|), when in fact they had been changed to
|Γ(x)| due to the change to __ieee754_gamma_r in the above patch.

Even today, wr_gamma.c says:

	double gamma_r(double x, int *signgamp) /* wrapper lgamma_r */

which is incorrect, as it calls __ieee754_gamma_r.

I see a couple of options here:

 1) Leave things as they are. This makes the 'gamma' family of functions
    incompatible with all other libc implementations and creates a trap
    for applications expecting either lgamma or tgamma values.

 2) Finish the work started in 2002 and make the 'gamma' family
    compatible with their 'tgamma' equivalents. This is what my patch
    does. This makes newlib compatible with 4.4 BSD libc, and leaves
    'gamma' returning the same value anytime that value is non-negative.

 3) Revert the 2002 change so that all of the 'gamma' family return the
    same value as their lgamma equivalents. That would make gamma_r
    useful again, as an alias for lgamma_r. This would make newlib
    compatible with glibc.

 4) Remove the 'gamma' family completely. Given that the meaning of
    'gamma' differs between 4.4 BSD and glibc, it would be safest to
    force applications to select between tgamma and lgamma.

I chose option 2 because that offers the best API compatibility -- the
gamma functions return the same value for most arguments (any positive
arguments, and half of negative arguments).

Any code using the gamma_r functions are likely expecting that to return
lgamma_r (as documented in several places). I can't see how we can leave
that function in the library as-is responsibly.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[C Howland]

> ------------------------------
> *From:* Newlib <newlib-bounces@sourceware.org> on behalf of Keith Packard
> via Newlib <newlib@sourceware.org>
> *Sent:* Wednesday, August 26, 2020 1:03 PM
> *To:* newlib@sourceware.org <newlib@sourceware.org>
> *Subject:* [PATCH 0/3] libm: Clean up gamma functions
>
> EXTERNAL EMAIL - This email originated from outside of CACI. Do not click
> any links or attachments unless you recognize and trust the sender.
>
>
>
>
>
> This  series  cleans   up  some  API  oddities   and  inaccuracies  in
> errno/exception reporting  from the  various gamma  related functions.
> With this series applied, newlib now  matches the POSIX spec (and glibc)
> for  tgamma/lgamma  functions. I've  left  the  'gamma' functions  BSD
> compatible (tgamma instead of lgamma) to avoid changing the API.
>
> [PATCH 1/3] libm: Fix sign value returned from __ieee754_lgamma*_r(-0)
>
>         Following the C spec, gamma(-0) should be -INFINITY rather than
>         +INFINITY, so the *signgamp value in lgamma should be -1
>         rather than +1 for this case.
>
What C spec?  Neither C99 nor C11 say so that I see for lgamma().  POSIX
makes the direct statement "If x is a non-positive integer, a pole error
shall occur and lgamma(), lgammaf(), and lgammal() shall return +HUGE_VAL,
+HUGE_VALF, and +HUGE_VALL, respectively."  0, regardless of sign, falls
under that.
https://pubs.opengroup.org/onlinepubs/9699919799/functions/lgamma.html
Additionally, the Linux lgamma man page, using GLIBC, says the same thing
as quoted from POSIX (not word for word, but the identical result).
(POSIX does require the +-0 +-INFINITY for tgamma().  Did these get crossed
up?)
Craig

>
> [PATCH 2/3] libm: Remove __ieee754_gamma_r variants
>
>         Back in 2002, a patch changed the various gamma functions so
>         that the returned abs(gamma) instead of the ln of gamma. This
>         left the _r variants in place, with their extra parameter to
>         return the sign of the resulting value. This doesn't make any
>         sense to me and looks like it was just a mistake. I've changed
>         the __ieee754_gamma functions to apply the sign value locally
>         and eliminated the gamma_r versions as those are no longer useful.
>
> [PATCH 3/3] libm: Adjust errno/exception values for gamma/lgamma
>
>         There's a pretty detailed discussion in the commit message
>         here, the short version is that tgamma and lgamma are
>         specified with different errno/exception behaviour for
>         one class of inputs (negative integers).
>

Re: Fw: [PATCH 3/3] libm: Adjust errno/exception values for gamma/lgamma

[C Howland]

> ------------------------------
> *From:* Newlib <newlib-bounces@sourceware.org> on behalf of Keith Packard
> via Newlib <newlib@sourceware.org>
> *Sent:* Wednesday, August 26, 2020 1:03 PM
> *To:* newlib@sourceware.org <newlib@sourceware.org>
> *Subject:* [PATCH 3/3] libm: Adjust errno/exception values for
> gamma/lgamma
>
>
>
>
> This makes all of gamma/lgamma functions match POSIX (and glibc) for
> errno and exception values for both gamma and lgamma functions.
>
> The big change is to support the exception/errno value differences for
> tgamma/lgamma for negative integer arguments. tgamma produces
> EDOM/FE_INVALID while lgamma produces ERANGE/FE_DIVBYZERO.
>
> This was done by splitting the lowest level __ieee754_lgamma*_r
> functions apart. The new lower-level ___ieee754_lgamma*_r functions
> can now produce exceptions which are correct for either lgamma or
> tgamma and use the value in the signgamp parameter to select which
> mode to operate in.
>
> All functions now return EDOM/FE_INVALID for -INFINITY, although POSIX
> doesn't specify this behavior for lgamma. It does match glibc at
> least.
>
I think you're looking at an old GLIBC.  From my RHEL7 lgamma() man  page:
"       If the result overflows, a range error occurs, and the functions
return
       HUGE_VAL, HUGE_VALF, or HUGE_VALL, respectively, with the correct
math‐
       ematical sign.
...
BUGS
       In glibc 2.9 and earlier, when a pole error occurs,  errno  is  set
 to
       EDOM;  instead of the POSIX-mandated ERANGE.  Since version 2.10,
glibc
       does the right thing."
POSIX has only pole and range errors listed for lgamma(), both of which
return ERANGE, so it is definitely improper to return EDOM.  It should be
range, so ERANGE/FE_OVERFLOW is what it ought to be.

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

Re: Fw: [PATCH 3/3] libm: Adjust errno/exception values for gamma/lgamma

[Brian Inglis]

On 2020-08-27 11:55, C Howland via Newlib wrote:
> On Wednesday, August 26, 2020 1:03 PM, Keith Packard wrote:
>> This makes all of gamma/lgamma functions match POSIX (and glibc) for
>> errno and exception values for both gamma and lgamma functions.
>>
>> The big change is to support the exception/errno value differences for
>> tgamma/lgamma for negative integer arguments. tgamma produces
>> EDOM/FE_INVALID while lgamma produces ERANGE/FE_DIVBYZERO.
>>
>> This was done by splitting the lowest level __ieee754_lgamma*_r
>> functions apart. The new lower-level ___ieee754_lgamma*_r functions
>> can now produce exceptions which are correct for either lgamma or
>> tgamma and use the value in the signgamp parameter to select which
>> mode to operate in.
>>
>> All functions now return EDOM/FE_INVALID for -INFINITY, although POSIX
>> doesn't specify this behavior for lgamma. It does match glibc at
>> least.

> I think you're looking at an old GLIBC.  From my RHEL7 lgamma() man  page:
> "       If the result overflows, a range error occurs, and the functions
> return
>        HUGE_VAL, HUGE_VALF, or HUGE_VALL, respectively, with the correct
> math‐
>        ematical sign.
> ...
> BUGS
>        In glibc 2.9 and earlier, when a pole error occurs,  errno  is  set
>  to
>        EDOM;  instead of the POSIX-mandated ERANGE.  Since version 2.10,
> glibc
>        does the right thing."
> POSIX has only pole and range errors listed for lgamma(), both of which
> return ERANGE, so it is definitely improper to return EDOM.  It should be
> range, so ERANGE/FE_OVERFLOW is what it ought to be.

RHEL7 uses an *old glibc* 2.17 from 2012, compared to Fedora Rawhide with latest
glibc release 2.32.

Check out the specs at the tops of the *current glibc* sources of the gamma
functions to be found nearby online under:

	https://sourceware.org/git/?p=glibc.git;a=tree;f=sysdeps/ieee754/dbl-64

-- 
Take care. Thanks, Brian Inglis, Calgary, Alberta, Canada

This email may be disturbing to some readers as it contains
too much technical detail. Reader discretion is advised.
[Data in IEC units and prefixes, physical quantities in SI.]

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 1719 bytes --]

C Howland via Newlib <newlib@sourceware.org> writes:

>>         Following the C spec, gamma(-0) should be -INFINITY rather than
>>         +INFINITY, so the *signgamp value in lgamma should be -1
>>         rather than +1 for this case.
>>
> What C spec?  Neither C99 nor C11 say so that I see for lgamma().  POSIX
> makes the direct statement "If x is a non-positive integer, a pole error
> shall occur and lgamma(), lgammaf(), and lgammal() shall return +HUGE_VAL,
> +HUGE_VALF, and +HUGE_VALL, respectively."  0, regardless of sign, falls
> under that.

Right, C99/C11 doesn't mention signgam at all, so POSIX is free to
define it as it likes. That spec only says that it's value is undefined
when the parameter is a negative integer.

C99 says that , tgamma should return -INFINITY for -0. As newlib makes
'gamma' an alias for tgamma, gamma should probably return the same
thing. To make that work with our implementation of tgamma, (which uses
lgamma and the 'signgam' value), that signgam value should be -1.

> https://pubs.opengroup.org/onlinepubs/9699919799/functions/lgamma.html
> Additionally, the Linux lgamma man page, using GLIBC, says the same thing
> as quoted from POSIX (not word for word, but the identical result).
> (POSIX does require the +-0 +-INFINITY for tgamma().  Did these get crossed
> up?)
> Craig

I think we're probably just confusing gamma between lgamma and tgamma --
newlib defines gamma 'oddly' (fixed in 2/3).

I think the only problematic patch in the series is in changing the
semantics of gamma/gammaf and removing the gamma_r/gammaf_r -- we need
to decide how to recover from what looks like a mistake introduced 18
years ago.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Brian Inglis]

On 2020-08-27 17:59, Keith Packard via Newlib wrote:
> C Howland via Newlib <newlib@sourceware.org> writes:
> 
>>> Following the C spec, gamma(-0) should be -INFINITY rather than 
>>> +INFINITY, so the *signgamp value in lgamma should be -1 rather than +1
>>> for this case.
>>>
>> What C spec?  Neither C99 nor C11 say so that I see for lgamma().  POSIX
>> makes the direct statement "If x is a non-positive integer, a pole error
>> shall occur and lgamma(), lgammaf(), and lgammal() shall return +HUGE_VAL,
>> +HUGE_VALF, and +HUGE_VALL, respectively."  0, regardless of sign, falls
>> under that.
> 
> Right, C99/C11 doesn't mention signgam at all, so POSIX is free to
> define it as it likes. That spec only says that it's value is undefined
> when the parameter is a negative integer.
> 
> C99 says that , tgamma should return -INFINITY for -0. As newlib makes
> 'gamma' an alias for tgamma, gamma should probably return the same
> thing. To make that work with our implementation of tgamma, (which uses
> lgamma and the 'signgam' value), that signgam value should be -1.

Last C17 public FDIS N2176 is available via:

https://web.archive.org/web/20181230041359if_/http://www.open-std.org/jtc1/sc22/wg14/www/abq/c17_updated_proposed_fdis.pdf

see 7.12.8.3 lgamma and 7.12.8.4 tgamma on pp181-182, 7.26 tgmath #5 lgamma,
tgamma pp272-273, 7.31.1 complex clgamma, ctgamma p332, F10.5.3 lgamma and
F.10.5.4 tgamma.

>> https://pubs.opengroup.org/onlinepubs/9699919799/functions/lgamma.html
>> Additionally, the Linux lgamma man page, using GLIBC, says the same thing
>> as quoted from POSIX (not word for word, but the identical result).
>> (POSIX does require the +-0 +-INFINITY for tgamma().  Did these get crossed
>> up?)
>> Craig
> 
> I think we're probably just confusing gamma between lgamma and tgamma --
> newlib defines gamma 'oddly' (fixed in 2/3).
> 
> I think the only problematic patch in the series is in changing the
> semantics of gamma/gammaf and removing the gamma_r/gammaf_r -- we need
> to decide how to recover from what looks like a mistake introduced 18
> years ago.

-- 
Take care. Thanks, Brian Inglis, Calgary, Alberta, Canada

This email may be disturbing to some readers as it contains
too much technical detail. Reader discretion is advised.
[Data in IEC units and prefixes, physical quantities in SI.]

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 5218 bytes --]

Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:

> Last C17 public FDIS N2176 is available via:
>
> https://web.archive.org/web/20181230041359if_/http://www.open-std.org/jtc1/sc22/wg14/www/abq/c17_updated_proposed_fdis.pdf
>
> see 7.12.8.3 lgamma and 7.12.8.4 tgamma on pp181-182, 7.26 tgmath #5 lgamma,
> tgamma pp272-273, 7.31.1 complex clgamma, ctgamma p332, F10.5.3 lgamma and
> F.10.5.4 tgamma.

Thanks for the link -- the difference for lgamma is that in C99:

        A range error occurs if x is too large. A range error may occur
        if x is a negative integer or zero.

While in C17:

        A range error occurs if positive x is too large. A pole error
        may occur if x is a negative integer or zero.

Frustratingly, tgamma still returns different errors. In C99:

        A domain error or range error may occur if x is a negative
        integer or zero. A range error may occur if the magnitude of x
        is too large or too small.

While in C17:

        A domain error or pole error may occur if x is a negative
        integer or zero. A range error occurs if the magnitude of x is
        too large and may occur if the magnitude of x is too small.

C17 appears to have added this new 'pole' error, which seems to match
what the IEEE 754 calls a 'divideByZero exception'. Let's see if my
understanding of the various errors is correct:

IEEE            C17             POSIX   errno   Exception
divideByZero    Pole            Pole    ERANGE  FE_DIVBYZERO
invalid         Domain          Domain  EDOM    FE_INVALID
overflow        Range           Range   ERANGE  FE_OVERFLOW

POSIX is more specific than C17, but it does appear to have been updated
to follow that spec, including the term 'Pole Error' which wasn't in
C99.

As my tests all refer to the POSIX spec (which defines the
errno/exception values), I think they're correct. Here's the set of
values I'm testing for both tgamma and lgamma (I'm not testing gamma
currently; we need to sort out what it should be doing first):

Value                 tgamma                               lgamma
         POSIX   result   errno  exception      POSIX   result   errno   exception

+0       Pole   +INFINITY ERANGE FE_DIVBYZERO   Pole   +INFINITY ERANGE FE_DIVBYZERO
-0       Pole   -INFINITY ERANGE FE_DIVBYZERO   Pole   +INFINITY ERANGE FE_DIVBYZERO
-1       Domain    NAN    EDOM   FE_INVALID     Pole   +INFINITY ERANGE FE_DIVBYZERO
+3e38f   Range  +INFINITY ERANGE FE_OVERFLOW    Range  +INFINITY ERANGE FE_OVERFLOW
+1.7e308 Range  +INFINITY ERANGE FE_OVERFLOW    Range  +INFINITY ERANGE FE_OVERFLOW
-3e38f   Domain    NAN    EDOM   FE_INVALID     Pole   +INFINITY ERANGE FE_DIVBYZERO
-1.7e308 Domain    NAN    EDOM   FE_INVALID     Pole   +INFINITY ERANGE FE_DIVBYZERO
+inf     None   +INFINITY 0      0              None   +INFINITY 0      0
-inf     Domain    NAN    EDOM   FE_INVALID     None   +INFINITY 0      0

I have rationalized the differences in this way:

        Domain error means there is no defined limit for the function
        approaching the parameter. Pole error means that there *is* a
        defined limit in this case.

        For non-positive inputs, the limit of the value of the tgamma
        function depends on whether that limit is approached from above
        or below.

        For zero, we can separate out 'approach from below' and
        'approach from above' with the sign -- -0 might well be the
        result of an underflow from below, hence returning the value of
        the function as approached from below makes "sense", in this
        case -INFINITY. Similarly, +0 might well be underflow from
        above, so the function can return +INFINITY. In both cases, as
        the limit is well defined, this represents a Pole error rather
        than a Domain error, hence the differences in errno and exceptions.
        
        For negative integers, as the limit depends on whether the value
        is approached from below or above and we have no way of
        distinguishing these two cases from the input, the limit is not
        well defined, so we have a Domain error and return
        NAN/EDOM/FE_INVALID.

        For the lgamma function, the limit for non-positive integers is
        always +INFINITY because it returns ln(|Γ(x)|). This makes it
        well defined and hence a Pole error rather than a Domain error

        Huge negative values (-1e308) are negative integers, and hence
        treated that way, generating a Domain error for tgamma and Range
        error for lgamma

        +INFINITY input and +INFINITY output generate *no* error
        presumably because there's an assumption that the computation
        which resulted in +INFINITY being provided to these functions
        already raised a FE_OVERFLOW exception, and perhaps an ERANGE
        errno.

Please review the table above to see if you agree about what errno and
exception values each of the provided inputs should generate. I'd really
like to know if you can suggest additional test cases to use to validate
the correctness of these functions.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Brian Inglis]

[-- Attachment #1.1: Type: text/plain, Size: 6370 bytes --]

On 2020-08-27 21:13, Keith Packard wrote:
> Brian Inglis writes:
> 
>> Last C17 public FDIS N2176 is available via:
>>
>> https://web.archive.org/web/20181230041359if_/http://www.open-std.org/jtc1/sc22/wg14/www/abq/c17_updated_proposed_fdis.pdf
>>
>> see 7.12.8.3 lgamma and 7.12.8.4 tgamma on pp181-182, 7.26 tgmath #5 lgamma,
>> tgamma pp272-273, 7.31.1 complex clgamma, ctgamma p332, F10.5.3 lgamma and
>> F.10.5.4 tgamma.
> 
> Thanks for the link -- the difference for lgamma is that in C99:
> 
>         A range error occurs if x is too large. A range error may occur
>         if x is a negative integer or zero.
> 
> While in C17:
> 
>         A range error occurs if positive x is too large. A pole error
>         may occur if x is a negative integer or zero.
> 
> Frustratingly, tgamma still returns different errors. In C99:
> 
>         A domain error or range error may occur if x is a negative
>         integer or zero. A range error may occur if the magnitude of x
>         is too large or too small.
> 
> While in C17:
> 
>         A domain error or pole error may occur if x is a negative
>         integer or zero. A range error occurs if the magnitude of x is
>         too large and may occur if the magnitude of x is too small.
> 
> C17 appears to have added this new 'pole' error, which seems to match
> what the IEEE 754 calls a 'divideByZero exception'. Let's see if my
> understanding of the various errors is correct:
> 
> IEEE            C17             POSIX   errno   Exception
> divideByZero    Pole            Pole    ERANGE  FE_DIVBYZERO
> invalid         Domain          Domain  EDOM    FE_INVALID
> overflow        Range           Range   ERANGE  FE_OVERFLOW
> 
> POSIX is more specific than C17, but it does appear to have been updated
> to follow that spec, including the term 'Pole Error' which wasn't in
> C99.
> 
> As my tests all refer to the POSIX spec (which defines the
> errno/exception values), I think they're correct. Here's the set of
> values I'm testing for both tgamma and lgamma (I'm not testing gamma
> currently; we need to sort out what it should be doing first):
> 
> Value                 tgamma                               lgamma
>          POSIX   result   errno  exception      POSIX   result   errno   exception
> 
> +0       Pole   +INFINITY ERANGE FE_DIVBYZERO   Pole   +INFINITY ERANGE FE_DIVBYZERO
> -0       Pole   -INFINITY ERANGE FE_DIVBYZERO   Pole   +INFINITY ERANGE FE_DIVBYZERO
> -1       Domain    NAN    EDOM   FE_INVALID     Pole   +INFINITY ERANGE FE_DIVBYZERO
> +3e38f   Range  +INFINITY ERANGE FE_OVERFLOW    Range  +INFINITY ERANGE FE_OVERFLOW
> +1.7e308 Range  +INFINITY ERANGE FE_OVERFLOW    Range  +INFINITY ERANGE FE_OVERFLOW
> -3e38f   Domain    NAN    EDOM   FE_INVALID     Pole   +INFINITY ERANGE FE_DIVBYZERO
> -1.7e308 Domain    NAN    EDOM   FE_INVALID     Pole   +INFINITY ERANGE FE_DIVBYZERO
> +inf     None   +INFINITY 0      0              None   +INFINITY 0      0
> -inf     Domain    NAN    EDOM   FE_INVALID     None   +INFINITY 0      0
> 
> I have rationalized the differences in this way:
> 
>         Domain error means there is no defined limit for the function
>         approaching the parameter. Pole error means that there *is* a
>         defined limit in this case.
> 
>         For non-positive inputs, the limit of the value of the tgamma
>         function depends on whether that limit is approached from above
>         or below.
> 
>         For zero, we can separate out 'approach from below' and
>         'approach from above' with the sign -- -0 might well be the
>         result of an underflow from below, hence returning the value of
>         the function as approached from below makes "sense", in this
>         case -INFINITY. Similarly, +0 might well be underflow from
>         above, so the function can return +INFINITY. In both cases, as
>         the limit is well defined, this represents a Pole error rather
>         than a Domain error, hence the differences in errno and exceptions.
>         
>         For negative integers, as the limit depends on whether the value
>         is approached from below or above and we have no way of
>         distinguishing these two cases from the input, the limit is not
>         well defined, so we have a Domain error and return
>         NAN/EDOM/FE_INVALID.
> 
>         For the lgamma function, the limit for non-positive integers is
>         always +INFINITY because it returns ln(|Γ(x)|). This makes it
>         well defined and hence a Pole error rather than a Domain error
> 
>         Huge negative values (-1e308) are negative integers, and hence
>         treated that way, generating a Domain error for tgamma and Range
>         error for lgamma
> 
>         +INFINITY input and +INFINITY output generate *no* error
>         presumably because there's an assumption that the computation
>         which resulted in +INFINITY being provided to these functions
>         already raised a FE_OVERFLOW exception, and perhaps an ERANGE
>         errno.
> 
> Please review the table above to see if you agree about what errno and
> exception values each of the provided inputs should generate. I'd really
> like to know if you can suggest additional test cases to use to validate
> the correctness of these functions.

I don't get so don't do maths above my comfort level: but it is possible that
ORs in the specs could mean that some important implementations (e.g BSD and
glibc) may differ in errors reported depending on which, possibly different,
exceptional input values are passed. You could check the BSD and glibc specs and
tests for these functions to see what the expected outputs are for what inputs,
what exceptional values generate which errno values, and whether all
possibilities are allowed by all the specs. POSIX always defers to C where the
latter defines library implementations, but as long as they do not contradict
normative content, POSIX may add additional conditions or requirements.

-- 
Take care. Thanks, Brian Inglis, Calgary, Alberta, Canada

This email may be disturbing to some readers as it contains
too much technical detail. Reader discretion is advised.
[Data in IEC units and prefixes, physical quantities in SI.]


[-- Attachment #2: OpenPGP digital signature --]
[-- Type: application/pgp-signature, Size: 228 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Corinna Vinschen]

On Aug 27 10:05, Keith Packard via Newlib wrote:
> Corinna Vinschen via Newlib <newlib@sourceware.org> writes:
> 
> > Nevertheless, the symbols have to be exported for backward compatibility.
> > So you're saying aliasing gamma_r to lgamma_r and gammaf_r to lgammaf_r
> > in the DLLs export table would be sufficient?
> 
> It depends if you want the pre-2002 functionality or the post-2002
> functionality.

The important point is that gamma_r and gammaf_r are exported by Cygwin
and must be kept available, otherwise building the DLL fails.  If they
returned the wrong vaslue, that's a bug and needs fixing.

Whether gamma_r and gammaf_r are still functions on their own, or if
they are just defined as aliases to lgamma_r and lgammaf_r in
winsup/cygwin/common.din, as in

  gamma_r = lgamma_r NOSIGFE
  gammaf_r = lgammaf_r NOSIGFE

is up to you.


Corinna

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Corinna Vinschen]

On Aug 28 10:19, Corinna Vinschen via Newlib wrote:
> On Aug 27 10:05, Keith Packard via Newlib wrote:
> > Corinna Vinschen via Newlib <newlib@sourceware.org> writes:
> > 
> > > Nevertheless, the symbols have to be exported for backward compatibility.
> > > So you're saying aliasing gamma_r to lgamma_r and gammaf_r to lgammaf_r
> > > in the DLLs export table would be sufficient?
> > 
> > It depends if you want the pre-2002 functionality or the post-2002
> > functionality.
> 
> The important point is that gamma_r and gammaf_r are exported by Cygwin
> and must be kept available, otherwise building the DLL fails.  If they
> returned the wrong vaslue, that's a bug and needs fixing.
> 
> Whether gamma_r and gammaf_r are still functions on their own, or if
> they are just defined as aliases to lgamma_r and lgammaf_r in
> winsup/cygwin/common.din, as in
> 
>   gamma_r = lgamma_r NOSIGFE
>   gammaf_r = lgammaf_r NOSIGFE
> 
> is up to you.

On second thought, given gamma_r and gammaf_r are BSD functions,
it might be a good idea to keep the wrapper functions for other
targets as well.

Any input from the RTEMS guys, perhaps?


Thanks,
Corinna

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 1366 bytes --]

Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:

> I don't get so don't do maths above my comfort level: but it is possible that
> ORs in the specs could mean that some important implementations (e.g BSD and
> glibc) may differ in errors reported depending on which, possibly different,
> exceptional input values are passed. You could check the BSD and glibc specs and
> tests for these functions to see what the expected outputs are for what inputs,
> what exceptional values generate which errno values, and whether all
> possibilities are allowed by all the specs. POSIX always defers to C where the
> latter defines library implementations, but as long as they do not contradict
> normative content, POSIX may add additional conditions or
> requirements.

Yeah, my testing could be more permissive to allow any valid
implementation to pass. Right now, I'm interested in making newlib
compliant with the relevant specs and, where that doesn't conflict,
compatible with glibc. And that means having hard requirements that
match glibc -- I haven't found any places where glibc gets the wrong
answer at least...

There's one place where that's a bit hard -- the gamma function has
changed meaning over time and the current newlib gamma implementation
is closer to the BSD one than glibc, but implements neither exactly.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Joseph Myers]

On Thu, 27 Aug 2020, Keith Packard via Newlib wrote:

>         +INFINITY input and +INFINITY output generate *no* error
>         presumably because there's an assumption that the computation
>         which resulted in +INFINITY being provided to these functions
>         already raised a FE_OVERFLOW exception, and perhaps an ERANGE
>         errno.

The basic principles for infinite results from IEEE floating-point 
operations and libm functions are that overflow is when infinity is an 
*inexact* result from rounding (or, in some rounding modes, the largest 
finite value results from rounding but rounding with infinite exponent 
range would have produced an exponent that could not be represented in the 
actual format).  Divide by zero is when infinity is an *exact* result from 
*finite* arguments.  If infinity is an exact result and at least one 
argument is infinite, there is no error or exception (unless it's some of 
the new minimum/maximum operations in IEEE 754-2019 and the other argument 
is a signaling NaN).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 800 bytes --]

Joseph Myers <joseph@codesourcery.com> writes:

> If infinity is an exact result and at least one 
> argument is infinite, there is no error or exception (unless it's some of 
> the new minimum/maximum operations in IEEE 754-2019 and the other argument 
> is a signaling NaN).

Thanks for the clarification. That all seems consistent, and I think I
can now (mostly) understand how the POSIX error values relate back to
the IEEE 754-2019 specification.

I do wonder if POSIX will add the new minimum/maximum operations; they
seem more like what I expect wrt NaN handling than the existing
minimumNumber/fmin maximumNumber/fmax operations.

Any other concerns about the testing values I'm using? Should we adopt
the newlib patches necessary to have it pass those tests?

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: Fw: [PATCH 0/3] libm: Clean up gamma functions

[Joseph Myers]

On Fri, 28 Aug 2020, Keith Packard via Newlib wrote:

> I do wonder if POSIX will add the new minimum/maximum operations; they
> seem more like what I expect wrt NaN handling than the existing
> minimumNumber/fmin maximumNumber/fmax operations.

C bindings for the new operations are given in 
http://www.open-std.org/jtc1/sc22/wg14/www/docs/n2532.pdf - but the target 
completion date for the next revision of POSIX (2022) is before that for 
C2x (publication aiming for 2023), so the next revision of POSIX is more 
likely to have only features from the 2011/2018 edition of ISO C.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Fabian Schriever]

Hi Keith,

We welcome your efforts to clean up - and correct the error return 
values of - the gamma/lgamma/tgamma families.

Regarding Corinna's concern with lgamma_r/gamma_r being BSD-functions:

> You can't do  that.  gamma_r/gammaf_r/lgamma_r/lgammaf_r are BSD
> functions. They have been exported by Cygwin since 2001.  The entry
> points need  to be kept available with unchanged semantics.

We would favor the removal of all non C/POSIX(+XSI)-standard functions 
from the interface.
The reason to regard lgamma_r and gamma_r as "BSD-functions" comes from 
the fact that BSD (same as newlib) took fdlibm as the base for its libm 
back in the early 90s when some non-standard functions were part of that 
library.
We would encourage the use of only the C/POSIX(+XSI)-standard functions 
as the only way to get rid of the confusing semantics of the historical 
function interfaces otherwise the problem will only perpetuate into the 
future.

If something was changed in 2002, that is working incorrectly and no one 
found out until now, that is also not part of any standard, it suggests 
that no one is actually using it and should be able to be safely 
removed. Does Newlib have a policy to remove elements?

An interesting discussion about the standard lgamma/tgamma functions 
would be to discuss accuracy improvements in line with the glibc 
improvements from Joseph Myers (see 
https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=050f29c18873ec05ba04a4034bed8cb3f6ae4463).

Best regards,
Fabian

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 2033 bytes --]

Fabian Schriever <fabian.schriever@gtd-gmbh.de> writes:

> Hi Keith,
>
> We welcome your efforts to clean up - and correct the error return 
> values of - the gamma/lgamma/tgamma families.

Thanks! I'm doing continuous integration testing as a part of embedded
toolchain support for RISC-V as part of my SiFive dayjob; fixing bugs
like this is part of that process.

> We would favor the removal of all non C/POSIX(+XSI)-standard functions 
> from the interface.

Oh, that's probably the best idea of all. Applications should not be
using 'gamma' at all given the different definitions over time and
space.

Any thoughts about the newlib __ieee754 interfaces? Those are
essentially the same as the C/POSIX interfaces but do not use errno or
other global variables, reporting exceptions only through the fenv API.

> If something was changed in 2002, that is working incorrectly and no one 
> found out until now, that is also not part of any standard, it suggests 
> that no one is actually using it and should be able to be safely 
> removed. Does Newlib have a policy to remove elements?

I don't think it's 'newlib' which would need any policy; newlib is used
downstream in a wide variety of projects, including cygwin and picolibc,
which may have separate policies. Cygwin has binary interface
definitions, changing those could affect applications there.

I'm using newlib as part of picolibc which is used for embedded
toolchains where removing things from the ABI to fix bugs would be just
fine.

> An interesting discussion about the standard lgamma/tgamma functions 
> would be to discuss accuracy improvements in line with the glibc 
> improvements from Joseph Myers (see 
> https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=050f29c18873ec05ba04a4034bed8cb3f6ae4463).

I read through that patch and agree that it would be nice to incorporate
something similar. We need to be cautious as code cannot be directly
brought in from glibc due to licensing differences.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Corinna Vinschen]

On Sep  1 10:23, Keith Packard via Newlib wrote:
> Fabian Schriever <fabian.schriever@gtd-gmbh.de> writes:
> 
> > Hi Keith,
> >
> > We welcome your efforts to clean up - and correct the error return 
> > values of - the gamma/lgamma/tgamma families.
> 
> Thanks! I'm doing continuous integration testing as a part of embedded
> toolchain support for RISC-V as part of my SiFive dayjob; fixing bugs
> like this is part of that process.
> 
> > We would favor the removal of all non C/POSIX(+XSI)-standard functions 
> > from the interface.
> 
> Oh, that's probably the best idea of all. Applications should not be
> using 'gamma' at all given the different definitions over time and
> space.
> 
> Any thoughts about the newlib __ieee754 interfaces? Those are
> essentially the same as the C/POSIX interfaces but do not use errno or
> other global variables, reporting exceptions only through the fenv API.
> 
> > If something was changed in 2002, that is working incorrectly and no one 
> > found out until now, that is also not part of any standard, it suggests 
> > that no one is actually using it and should be able to be safely 
> > removed. Does Newlib have a policy to remove elements?
> 
> I don't think it's 'newlib' which would need any policy; newlib is used
> downstream in a wide variety of projects, including cygwin and picolibc,
> which may have separate policies. Cygwin has binary interface
> definitions, changing those could affect applications there.

Removing interfaces is not an option for Cygwin.  If you remove
functions from newlib exported as symbols in Cygwin, you must provide
replacement interfaces for Cygwin alone, otherwise Cygwin won't build
anymore.  And one step further, removing interfaces from the list of
exported symbols of the Cygwin DLL will break user space and that's
simply a no-no.


Corinna

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 824 bytes --]

Corinna Vinschen via Newlib <newlib@sourceware.org> writes:

> Removing interfaces is not an option for Cygwin.  If you remove
> functions from newlib exported as symbols in Cygwin, you must provide
> replacement interfaces for Cygwin alone, otherwise Cygwin won't build
> anymore.  And one step further, removing interfaces from the list of
> exported symbols of the Cygwin DLL will break user space and that's
> simply a no-no.

Then tell us what you want gamma and gamma_r to do for Cygwin. They're
non-standard interfaces for POSIX, C99 and C17, so you get to
pick.

Removing them for non-Cygwin uses appears to be the safest choice as
that will (intentionally) break applications which were expecting them
to implement either the BSD or Linux behavior, when in fact they do
neither (!).

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Corinna Vinschen]

On Sep  2 13:37, Keith Packard via Newlib wrote:
> Corinna Vinschen via Newlib <newlib@sourceware.org> writes:
> 
> > Removing interfaces is not an option for Cygwin.  If you remove
> > functions from newlib exported as symbols in Cygwin, you must provide
> > replacement interfaces for Cygwin alone, otherwise Cygwin won't build
> > anymore.  And one step further, removing interfaces from the list of
> > exported symbols of the Cygwin DLL will break user space and that's
> > simply a no-no.
> 
> Then tell us what you want gamma and gamma_r to do for Cygwin.

I did: https://sourceware.org/pipermail/newlib/2020/017946.html


Thanks,
Corinna

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Brian Inglis]

On 2020-09-03 02:04, Corinna Vinschen via Newlib wrote:
> On Sep  2 13:37, Keith Packard via Newlib wrote:
>> Corinna Vinschen via Newlib <newlib@sourceware.org> writes:
>>
>>> Removing interfaces is not an option for Cygwin.  If you remove
>>> functions from newlib exported as symbols in Cygwin, you must provide
>>> replacement interfaces for Cygwin alone, otherwise Cygwin won't build
>>> anymore.  And one step further, removing interfaces from the list of
>>> exported symbols of the Cygwin DLL will break user space and that's
>>> simply a no-no.
>>
>> Then tell us what you want gamma and gamma_r to do for Cygwin.
> 
> I did: https://sourceware.org/pipermail/newlib/2020/017946.html

FYI docs/spec:
https://sca.uwaterloo.ca/coldfire/gcc-doc/docs/libm_21.html#SEC21
http://www.ece.ualberta.ca/~cmpe401/docs/coldfire/libm.pdf#page=20
https://ftp.rtems.org/pub/rtems/docs/3.2.0/libm-3.2.0.ps - see p.18

-- 
Take care. Thanks, Brian Inglis, Calgary, Alberta, Canada

This email may be disturbing to some readers as it contains
too much technical detail. Reader discretion is advised.
[Data in IEC units and prefixes, physical quantities in SI.]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 950 bytes --]

Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:

>>> Then tell us what you want gamma and gamma_r to do for Cygwin.
>> 
>> I did: https://sourceware.org/pipermail/newlib/2020/017946.html

Oops! I appear to have missed that mail.

> FYI docs/spec:
> https://sca.uwaterloo.ca/coldfire/gcc-doc/docs/libm_21.html#SEC21
> http://www.ece.ualberta.ca/~cmpe401/docs/coldfire/libm.pdf#page=20
> https://ftp.rtems.org/pub/rtems/docs/3.2.0/libm-3.2.0.ps - see p.18

That's quite helpful actually and says that cygwin's gamma matches
glibc, which makes the change from 2002 just a bug. Switching gamma back
to lgamma is quite simple.

Here's what I did:

 1. Remove the gamma* function implementations and make them explicit
    aliases to lgamma*.

 2. Change the name of the __ieee754_gamma* functions to
    __ieee754_tgamma* and remove the _ieee754_gamma*_r variants

I've sent a patch doing this to the list.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Brian Inglis]

[-- Attachment #1.1: Type: text/plain, Size: 2109 bytes --]

On 2020-09-03 15:25, Keith Packard wrote:
> Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:
> 
>>>> Then tell us what you want gamma and gamma_r to do for Cygwin.
>>>
>>> I did: https://sourceware.org/pipermail/newlib/2020/017946.html
> 
> Oops! I appear to have missed that mail.
> 
>> FYI docs/spec:
>> https://sca.uwaterloo.ca/coldfire/gcc-doc/docs/libm_21.html#SEC21
>> http://www.ece.ualberta.ca/~cmpe401/docs/coldfire/libm.pdf#page=20
>> https://ftp.rtems.org/pub/rtems/docs/3.2.0/libm-3.2.0.ps - see p.18
> 
> That's quite helpful actually and says that cygwin's gamma matches
> glibc, which makes the change from 2002 just a bug. Switching gamma back
> to lgamma is quite simple.
> 
> Here's what I did:
> 
>  1. Remove the gamma* function implementations and make them explicit
>     aliases to lgamma*.
> 
>  2. Change the name of the __ieee754_gamma* functions to
>     __ieee754_tgamma* and remove the _ieee754_gamma*_r variants
> 
> I've sent a patch doing this to the list.

I don't think you can "remove the _ieee754_gamma*_r variants" as I believe
Cygwin has to keep these:

$ nm --defined-only --extern-only /bin/cygwin1.dll | fgrep gamma
00000001801a76c0 T __ieee754_gamma_r
00000001801a7f90 T __ieee754_gammaf_r
00000001801a76e0 T __ieee754_lgamma_r
00000001801a7fb0 T __ieee754_lgammaf_r
000000018018a8c0 T __lgammal_r
000000018018ce10 T __tgammal_r
000000018019f690 T gamma
00000001801a0a40 T gamma_r
00000001801a0200 T gammaf
00000001801a0b80 T gammaf_r
000000018019fb00 T lgamma
00000001801a0ae0 T lgamma_r
00000001801a0660 T lgammaf
00000001801a0c20 T lgammaf_r
000000018018ae10 T lgammal
000000018018ae90 T lgammal_r
000000018019fe50 T tgamma
00000001801a0970 T tgammaf
000000018018d290 T tgammal

although I don't have the Windows utilities to check the exports directly.

-- 
Take care. Thanks, Brian Inglis, Calgary, Alberta, Canada

This email may be disturbing to some readers as it contains
too much technical detail. Reader discretion is advised.
[Data in IEC units and prefixes, physical quantities in SI.]


[-- Attachment #2: OpenPGP digital signature --]
[-- Type: application/pgp-signature, Size: 228 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 344 bytes --]

Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:

> I don't think you can "remove the _ieee754_gamma*_r variants" as I believe
> Cygwin has to keep these:

I added a note to that effect in my patch email; if someone can let me
know, I'll add them back in. I'd prefer to leave them out if they aren't
needed though.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Brian Inglis]

[-- Attachment #1.1: Type: text/plain, Size: 995 bytes --]

On 2020-09-03 18:01, Keith Packard wrote:
> Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:
> 
>> I don't think you can "remove the _ieee754_gamma*_r variants" as I believe
>> Cygwin has to keep these:
> 
> I added a note to that effect in my patch email; if someone can let me
> know, I'll add them back in. I'd prefer to leave them out if they aren't
> needed though.

$ fgrep gamma newlib-cygwin/winsup/cygwin/common.din

and

$ fgrep gamma newlib-cygwin/x86_64-pc-cygwin/winsup/cygwin/cygwin.def

have the same 13 function names:

gamma
gamma_r
gammaf
gammaf_r
lgamma
lgamma_r
lgammaf
lgammaf_r
lgammal
lgammal_r
tgamma
tgammaf
tgammal

so these may be the official ABI, and the others don't matter?

-- 
Take care. Thanks, Brian Inglis, Calgary, Alberta, Canada

This email may be disturbing to some readers as it contains
too much technical detail. Reader discretion is advised.
[Data in IEC units and prefixes, physical quantities in SI.]


[-- Attachment #2: OpenPGP digital signature --]
[-- Type: application/pgp-signature, Size: 228 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 485 bytes --]

Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:

> gamma
> gamma_r
> gammaf
> gammaf_r
> lgamma
> lgamma_r
> lgammaf
> lgammaf_r
> lgammal
> lgammal_r
> tgamma
> tgammaf
> tgammal
>
> so these may be the official ABI, and the others don't matter?

I would hope so -- the __ieee754 functions seem mostly like an internal
implementation detail, but it's hard to know.

If this is true, then I think the patch as posted should be what we want.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Corinna Vinschen]

On Sep  3 18:37, Keith Packard via Newlib wrote:
> Brian Inglis <Brian.Inglis@SystematicSw.ab.ca> writes:
> 
> > gamma
> > gamma_r
> > gammaf
> > gammaf_r
> > lgamma
> > lgamma_r
> > lgammaf
> > lgammaf_r
> > lgammal
> > lgammal_r
> > tgamma
> > tgammaf
> > tgammal
> >
> > so these may be the official ABI, and the others don't matter?
> 
> I would hope so -- the __ieee754 functions seem mostly like an internal
> implementation detail, but it's hard to know.
> 
> If this is true, then I think the patch as posted should be what we want.

The ieee functions can go away.  The symbols from cygwin1.dll don't
count, the ones exported from /usr/lib/libcygwin.a do.  These are
defined in winsup/cygwin/{common.din,i686.din,x86_64.din}


Corinna

Re: [PATCH 2/3] libm: Remove __ieee754_gamma_r variants

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 317 bytes --]

Corinna Vinschen via Newlib <newlib@sourceware.org> writes:

> The ieee functions can go away.  The symbols from cygwin1.dll don't
> count, the ones exported from /usr/lib/libcygwin.a do.  These are
> defined in winsup/cygwin/{common.din,i686.din,x86_64.din}

Which matches what we had guessed.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]