Update gamma functions from code in picolibc

- fixes issue with inf sign when x is -0

diff --git a/newlib/libm/math/er_lgamma.c b/newlib/libm/math/er_lgamma.c
index 386a8a73bd585557db27b3ef735557b7d66e417d..5c88548fb015defaa35a4eeadcf0c2513966171f 100644 (file)
--- 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 3c6ba02afbbb9b44cea2d16b66e3982c3e7040a1..f88f6309264f6845902facaae9210104d5d50930 100644 (file)
--- 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 c090925109ca7c089b7e7d267aa1fd374ab4e94b..0f90dd4c6c03ef22c81da91c55c0c6bda0f68e08 100644 (file)
--- 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 88567b2eb5daa3f48c237991ebf6fb83e9101bf6..80aacf75733b15c002e79cb0cc8ee6bfcd961f06 100644 (file)
--- 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