Update gamma functions from code in picolibc

- fixes issue with inf sign when x is -0
This commit is contained in:
Jeff Johnston 2020-12-17 15:58:49 -05:00
parent 1dd3f69db5
commit b2f3d593ff
4 changed files with 47 additions and 28 deletions

View File

@ -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
@ -213,7 +213,8 @@ static double zero= 0.00000000000000000000e+00;
double __ieee754_lgamma_r(double x, int *signgamp)
#else
double __ieee754_lgamma_r(x, signgamp)
double x; int *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;
}

View File

@ -148,7 +148,8 @@ static float zero= 0.0000000000e+00;
float __ieee754_lgammaf_r(float x, int *signgamp)
#else
float __ieee754_lgammaf_r(x, signgamp)
float x; int *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;

View File

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

View File

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