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

[Keith Packard <keithp@keithp.com>]

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