public inbox for newlib-cvs@sourceware.org
help / color / mirror / Atom feed
* [newlib-cygwin] Update gamma functions from code in picolibc
@ 2020-12-17 21:23 Jeff Johnston
0 siblings, 0 replies; only message in thread
From: Jeff Johnston @ 2020-12-17 21:23 UTC (permalink / raw)
To: newlib-cvs
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
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2020-12-17 21:23 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2020-12-17 21:23 [newlib-cygwin] Update gamma functions from code in picolibc Jeff Johnston
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).