public inbox for newlib-cvs@sourceware.org help / color / mirror / Atom feed
From: Jeff Johnston <jjohnstn@sourceware.org> To: newlib-cvs@sourceware.org Subject: [newlib-cygwin] Update gamma functions from code in picolibc Date: Thu, 17 Dec 2020 21:23:51 +0000 (GMT) [thread overview] Message-ID: <20201217212351.88DEA3896815@sourceware.org> (raw) https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;h=b2f3d593ff6074b945eb0ff925c2740e98a88ae0 commit b2f3d593ff6074b945eb0ff925c2740e98a88ae0 Author: Jeff Johnston <jjohnstn@redhat.com> Date: Thu Dec 17 15:58:49 2020 -0500 Update gamma functions from code in picolibc - fixes issue with inf sign when x is -0 Diff: --- 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
reply other threads:[~2020-12-17 21:23 UTC|newest] Thread overview: [no followups] expand[flat|nested] mbox.gz Atom feed
Reply instructions: You may reply publicly to this message via plain-text email using any one of the following methods: * Save the following mbox file, import it into your mail client, and reply-to-all from there: mbox Avoid top-posting and favor interleaved quoting: https://en.wikipedia.org/wiki/Posting_style#Interleaved_style * Reply using the --to, --cc, and --in-reply-to switches of git-send-email(1): git send-email \ --in-reply-to=20201217212351.88DEA3896815@sourceware.org \ --to=jjohnstn@sourceware.org \ --cc=newlib-cvs@sourceware.org \ /path/to/YOUR_REPLY https://kernel.org/pub/software/scm/git/docs/git-send-email.html * If your mail client supports setting the In-Reply-To header via mailto: links, try the mailto: linkBe sure your reply has a Subject: header at the top and a blank line before the message body.
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).