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

Message ID 20200826170357.2551683-4-keithp@keithp.com
State New
Headers show
Series
  • libm: Clean up gamma functions
Related show

Commit Message

Torbjorn SVENSSON via Newlib Aug. 26, 2020, 5:03 p.m.
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(-)

-- 
2.28.0

Comments

Torbjorn SVENSSON via Newlib Aug. 27, 2020, 5:55 p.m. | #1
> ------------------------------

> *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>

>

>
Brian Inglis Aug. 27, 2020, 7:28 p.m. | #2
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.]

Patch

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 */