From b2f3d593ff6074b945eb0ff925c2740e98a88ae0 Mon Sep 17 00:00:00 2001
From: Jeff Johnston <jjohnstn@redhat.com>
Date: Thu, 17 Dec 2020 15:58:49 -0500
Subject: [PATCH] Update gamma functions from code in picolibc

- fixes issue with inf sign when x is -0
---
 newlib/libm/math/er_lgamma.c  | 29 ++++++++++++++++++++---------
 newlib/libm/math/erf_lgamma.c | 25 ++++++++++++++++++-------
 newlib/libm/math/w_tgamma.c   | 10 +++++-----
 newlib/libm/math/wf_tgamma.c  | 11 ++++-------
 4 files changed, 47 insertions(+), 28 deletions(-)

diff --git a/newlib/libm/math/er_lgamma.c b/newlib/libm/math/er_lgamma.c
index 386a8a73b..5c88548fb 100644
--- a/newlib/libm/math/er_lgamma.c
+++ b/newlib/libm/math/er_lgamma.c
@@ -12,9 +12,9 @@
  *
  */
 
-/* __ieee754_lgamma_r(x, signgamp)
+/* __ieee754_lgamma_r(x)
  * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
+ * with signgam for the sign of Gamma(x). 
  *
  * Method:
  *   1. Argument Reduction for 0 < x <= 8
@@ -212,8 +212,9 @@ static double zero=  0.00000000000000000000e+00;
 #ifdef __STDC__
 	double __ieee754_lgamma_r(double x, int *signgamp)
 #else
-	double __ieee754_lgamma_r(x,signgamp)
-	double x; int *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;
@@ -224,8 +225,14 @@ static double zero=  0.00000000000000000000e+00;
     /* purge off +-inf, NaN, +-0, and negative arguments */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
-	if(ix>=0x7ff00000) return x*x;
-	if((ix|lx)==0) return one/zero;
+	if(ix>=0x7ff00000) {
+	    return x*x;
+	}
+	if((ix|lx)==0) {
+	    if(hx<0)
+	        *signgamp = -1;
+	    return one/(x-x);
+	}
 	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
@@ -233,10 +240,13 @@ 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 */
+		return one/(x-x); /* -integer */
+	    }
 	    t = sin_pi(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) {
+		return one/(x-x); /* -integer */
+	    }
 	    nadj = __ieee754_log(pi/fabs(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
@@ -307,3 +317,4 @@ static double zero=  0.00000000000000000000e+00;
 	if(hx<0) r = nadj - r;
 	return r;
 }
+
diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
index 3c6ba02af..f88f63092 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -147,8 +147,9 @@ static float zero=  0.0000000000e+00;
 #ifdef __STDC__
 	float __ieee754_lgammaf_r(float x, int *signgamp)
 #else
-	float __ieee754_lgammaf_r(x,signgamp)
-	float x; int *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;
@@ -159,8 +160,14 @@ static float zero=  0.0000000000e+00;
     /* purge off +-inf, NaN, +-0, and negative arguments */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) return x*x;
-	if(ix==0) return one/zero;
+	if(ix>=0x7f800000) {
+	    return x*x;
+	}
+	if(ix==0) {
+	    if(hx<0)
+	        *signgamp = -1;
+	    return one/(x-x);
+	}
 	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
@@ -168,10 +175,14 @@ static float zero=  0.0000000000e+00;
 	    } else return -__ieee754_logf(x);
 	}
 	if(hx<0) {
-	    if(ix>=0x4b000000) 	/* |x|>=2**23, must be -integer */
-		return one/zero;
+	    if(ix>=0x4b000000) { 	/* |x|>=2**23, must be -integer */
+		return one/(x-x);
+	    }
 	    t = sin_pif(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) {
+		/* tgamma wants NaN instead of INFINITY */
+		return one/(x-x); /* -integer */
+	    }
 	    nadj = __ieee754_logf(pi/fabsf(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
diff --git a/newlib/libm/math/w_tgamma.c b/newlib/libm/math/w_tgamma.c
index c09092510..0f90dd4c6 100644
--- a/newlib/libm/math/w_tgamma.c
+++ b/newlib/libm/math/w_tgamma.c
@@ -33,11 +33,11 @@
 #else
 	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 */
-	  else
-	    return __kernel_standard(x,x,40); /* tgamma overflow */
+	if(!finite(y)) {
+	  if(x < 0.0 && floor(x)==x)
+	    errno = EDOM;
+	  else if (finite(x))
+	    errno = ERANGE;
 	}
 	return y;
 #endif
diff --git a/newlib/libm/math/wf_tgamma.c b/newlib/libm/math/wf_tgamma.c
index 88567b2eb..80aacf757 100644
--- a/newlib/libm/math/wf_tgamma.c
+++ b/newlib/libm/math/wf_tgamma.c
@@ -30,13 +30,10 @@
 #else
 	if(_LIB_VERSION == _IEEE_) return y;
 
-	if(!finitef(y)&&finitef(x)) {
-	  if(floorf(x)==x&&x<=(float)0.0)
-	    /* tgammaf pole */
-	    return (float)__kernel_standard((double)x,(double)x,141);
-	  else
-	    /* tgammaf overflow */
-	    return (float)__kernel_standard((double)x,(double)x,140);
+	if(x < 0.0 && floor(x)==x)
+	    errno = EDOM;
+	  else if (finite(x))
+	    errno = ERANGE;
 	}
 	return y;
 #endif