regression in Newlib 4.2.0 (lgammaf)

regression in Newlib 4.2.0 (lgammaf)

[Paul Zimmermann]

       Hi,

a third regression in Newlib 4.2.0 for single precision:

Checking lgamma with newlib-4.2.0.20211231
MPFR library: 4.1.0       
MPFR header:  4.1.0 (based on 4.1.0)
Checking function mylgammaf with MPFR_RNDN
libm wrong by up to 8.93e+43 ulp(s) for x=-0x1p-70
mylgamma      gives inf
mpfr_mylgamma gives 0x1.842994p+5

The largest error for lgammaf in Newlib 4.1.0 was 7.50e6 ulps
(for x=-0x1.3a7fcap+1).

Paul

[PATCH 1/2] Improve lgammaf range for very small cases

[Andoni Arregi]

The original cut for small arguments at |x|<2**-70 (copied from the
double version) produces that when computing nadj we get a subnormal
number for t*x and thus, the division of pi/subnormal will be INF and
the logarithm of it too, which is wrong as a result for lgammaf in this
range.
The proposed new limit seems to be safe and has been tested to
produce accurate results.
(Courtesy of Andreas Jung, ESA)
---
 newlib/libm/math/erf_lgamma.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/newlib/libm/math/erf_lgamma.c
b/newlib/libm/math/erf_lgamma.c
index f88f63092..84d02159b 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
 	       *signgamp = -1;
 	   return one/(x-x);
 	}
-	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
+	if(ix<0x30800000) {	/* |x|<2**-30, return -log(|x|) */
 	   if(hx<0) {
 	       *signgamp = -1;
 	       return -__ieee754_logf(-x);
-- 
2.35.1

Re: [PATCH 1/2] Improve lgammaf range for very small cases

[Corinna Vinschen]

Hi Andoni,

On Feb  9 17:28, Andoni Arregi wrote:
> The original cut for small arguments at |x|<2**-70 (copied from the
> double version) produces that when computing nadj we get a subnormal
> number for t*x and thus, the division of pi/subnormal will be INF and
> the logarithm of it too, which is wrong as a result for lgammaf in this
> range.
> The proposed new limit seems to be safe and has been tested to
> produce accurate results.
> (Courtesy of Andreas Jung, ESA)
> ---
>  newlib/libm/math/erf_lgamma.c | 2 +-
>  1 file changed, 1 insertion(+), 1 deletion(-)
> 
> diff --git a/newlib/libm/math/erf_lgamma.c
> b/newlib/libm/math/erf_lgamma.c

This is broken.  It should be on one line, so apparently your
mailer broke it via autowrapping lines.

> index f88f63092..84d02159b 100644
> --- a/newlib/libm/math/erf_lgamma.c
> +++ b/newlib/libm/math/erf_lgamma.c
> @@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
>  	       *signgamp = -1;
>  	   return one/(x-x);
>  	}
> -	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
> +	if(ix<0x30800000) {	/* |x|<2**-30, return -log(|x|) */
>  	   if(hx<0) {
>  	       *signgamp = -1;
>  	       return -__ieee754_logf(-x);
> -- 
> 2.35.1

But even if I fix this manually, your patch doesn't apply, neither with
`git am', nor with `patch -p1'.  I have to admit that I don't understand
why...

Another question, would you mind to send patch series with a cover letter,
please?


Thanks,
Corinna

[PATCH 0/2] Fix lgammaf and lgamma for small ranges and in [2, 3[

[Andoni Arregi]

lgammaf shows a wrong behavior for very small cases and also for cases
in the range [2.0, 3.0), this last problem also being present in
lgamma. These wrong behaviors also affect the tgammaf and tgamma
functions respectively. The provided patches fix these errors.

Andoni Arregi (2):
  Improve lgammaf range for very small cases
  Add a missing default case in lgamma

 newlib/libm/math/er_lgamma.c  | 1 +
 newlib/libm/math/erf_lgamma.c | 3 ++-
 2 files changed, 3 insertions(+), 1 deletion(-)

-- 
2.35.1

[PATCH 1/2] Improve lgammaf range for very small cases

[Andoni Arregi]

The original cut for small arguments at |x|<2**-70 (copied from the
double version) produces that when computing nadj we get a subnormal
number for t*x and thus, the division of pi/subnormal will be INF and
the logarithm of it too, which is wrong as a result for lgammaf in this
range.
The proposed new limit seems to be safe and has been tested to
produce accurate results.
(Courstesy of Andreas Jung, ESA)
---
 newlib/libm/math/erf_lgamma.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
index f88f63092..84d02159b 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
 	       *signgamp = -1;
 	   return one/(x-x);
 	}
-	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
+	if(ix<0x30800000) {	/* |x|<2**-30, return -log(|x|) */
 	   if(hx<0) {
 	       *signgamp = -1;
 	       return -__ieee754_logf(-x);
-- 
2.35.1

Re: [PATCH 1/2] Improve lgammaf range for very small cases

[Paul Zimmermann]

       Dear Andoni,

I have tested this patch: the maximal error for all binary32 inputs to lgammaf
decreases to 7.50e+06 ulps, which is still huge, but less huge than 8.93e+43...

Using RedHat newlib
MPFR library: 4.1.0       
MPFR header:  4.1.0 (based on 4.1.0)
Checking function mylgammaf with MPFR_RNDN
libm wrong by up to 7.50e+06 ulp(s) [7497618] for x=-0x1.3a7fcap+1
mylgamma      gives -0x1p-24
mpfr_mylgamma gives -0x1.e4cf24p-24
Total: errors=509423944 (11.91%) errors2=11684280 maxerr=7.50e+06 ulp(s)

Best regards,
Paul

> From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> Date: Thu, 10 Feb 2022 17:10:35 +0100
> Organization: GTD GmbH
> User-Agent: Evolution 3.42.3 
> 
> The original cut for small arguments at |x|<2**-70 (copied from the
> double version) produces that when computing nadj we get a subnormal
> number for t*x and thus, the division of pi/subnormal will be INF and
> the logarithm of it too, which is wrong as a result for lgammaf in this
> range.
> The proposed new limit seems to be safe and has been tested to
> produce accurate results.
> (Courstesy of Andreas Jung, ESA)
> ---
>  newlib/libm/math/erf_lgamma.c | 2 +-
>  1 file changed, 1 insertion(+), 1 deletion(-)
> 
> diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
> index f88f63092..84d02159b 100644
> --- a/newlib/libm/math/erf_lgamma.c
> +++ b/newlib/libm/math/erf_lgamma.c
> @@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
>  	       *signgamp = -1;
>  	   return one/(x-x);
>  	}
> -	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
> +	if(ix<0x30800000) {	/* |x|<2**-30, return -log(|x|) */
>  	   if(hx<0) {
>  	       *signgamp = -1;
>  	       return -__ieee754_logf(-x);
> -- 
> 2.35.1
> 
> 
> 

[PATCH 2/2] Add a missing default case in lgamma

[Andoni Arregi]

The missing default case leads to large errors for |x| in range
[2.0, 3.0[.
In case of i=2, i.e., in the range between [2,3[, the computation
"r += logf(z)" is not performed and hence the result was wrong.
The added default case resolves this issue.
(Courtesy of Andreas Jung, ESA)
---
 newlib/libm/math/er_lgamma.c  | 1 +
 newlib/libm/math/erf_lgamma.c | 1 +
 2 files changed, 2 insertions(+)

diff --git a/newlib/libm/math/er_lgamma.c b/newlib/libm/math/er_lgamma.c
index 5c88548fb..65727c6ab 100644
--- a/newlib/libm/math/er_lgamma.c
+++ b/newlib/libm/math/er_lgamma.c
@@ -302,6 +302,7 @@ static double zero=  0.00000000000000000000e+00;
 	   case 5: z *= (y+4.0);	/* FALLTHRU */
 	   case 4: z *= (y+3.0);	/* FALLTHRU */
 	   case 3: z *= (y+2.0);	/* FALLTHRU */
+      default:
 		   r += __ieee754_log(z); break;
 	   }
     /* 8.0 <= x < 2**58 */
diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
index 84d02159b..e7311cacf 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -238,6 +238,7 @@ static float zero=  0.0000000000e+00;
 	   case 5: z *= (y+(float)4.0);	/* FALLTHRU */
 	   case 4: z *= (y+(float)3.0);	/* FALLTHRU */
 	   case 3: z *= (y+(float)2.0);	/* FALLTHRU */
+      default:
 		   r += __ieee754_logf(z); break;
 	   }
     /* 8.0 <= x < 2**58 */
-- 
2.35.1

Re: [PATCH 2/2] Add a missing default case in lgamma

[Paul Zimmermann]

       Dear Andoni,

unless I did a mistake, I still get the same error (7.50e+06 ulps)
than before this patch (after applying the other one), when I do an
exhaustive search on all binary32 inputs:

Using RedHat newlib
MPFR library: 4.1.0       
MPFR header:  4.1.0 (based on 4.1.0)
Checking function mylgammaf with MPFR_RNDN
libm wrong by up to 7.50e+06 ulp(s) [7497618] for x=-0x1.3a7fcap+1
mylgamma      gives -0x1p-24
mpfr_mylgamma gives -0x1.e4cf24p-24
Total: errors=509423944 (11.91%) errors2=11684280 maxerr=7.50e+06 ulp(s)

What result do you get for x=-0x1.3a7fcap+1?

Best regards,
Paul

> From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> Date: Thu, 10 Feb 2022 17:12:40 +0100
> 
> The missing default case leads to large errors for |x| in range
> [2.0, 3.0[.
> In case of i=2, i.e., in the range between [2,3[, the computation
> "r += logf(z)" is not performed and hence the result was wrong.
> The added default case resolves this issue.
> (Courtesy of Andreas Jung, ESA)
> ---
>  newlib/libm/math/er_lgamma.c  | 1 +
>  newlib/libm/math/erf_lgamma.c | 1 +
>  2 files changed, 2 insertions(+)
> 
> diff --git a/newlib/libm/math/er_lgamma.c b/newlib/libm/math/er_lgamma.c
> index 5c88548fb..65727c6ab 100644
> --- a/newlib/libm/math/er_lgamma.c
> +++ b/newlib/libm/math/er_lgamma.c
> @@ -302,6 +302,7 @@ static double zero=  0.00000000000000000000e+00;
>  	   case 5: z *= (y+4.0);	/* FALLTHRU */
>  	   case 4: z *= (y+3.0);	/* FALLTHRU */
>  	   case 3: z *= (y+2.0);	/* FALLTHRU */
> +      default:
>  		   r += __ieee754_log(z); break;
>  	   }
>      /* 8.0 <= x < 2**58 */
> diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
> index 84d02159b..e7311cacf 100644
> --- a/newlib/libm/math/erf_lgamma.c
> +++ b/newlib/libm/math/erf_lgamma.c
> @@ -238,6 +238,7 @@ static float zero=  0.0000000000e+00;
>  	   case 5: z *= (y+(float)4.0);	/* FALLTHRU */
>  	   case 4: z *= (y+(float)3.0);	/* FALLTHRU */
>  	   case 3: z *= (y+(float)2.0);	/* FALLTHRU */
> +      default:
>  		   r += __ieee754_logf(z); break;
>  	   }
>      /* 8.0 <= x < 2**58 */
> -- 
> 2.35.1

Re: [PATCH 2/2] Add a missing default case in lgamma

[Andoni Arregi]

On Thu, 2022-02-10 at 18:29 +0100, Paul Zimmermann wrote:
>        Dear Andoni,
> 
> unless I did a mistake, I still get the same error (7.50e+06 ulps)
> than before this patch (after applying the other one), when I do an
> exhaustive search on all binary32 inputs:
> 
You didn't make any mistake. I did so when porting back our fixes.
As you say, the first patch does improve the lgammaf function from about
8.93e+43 ULP to 7.50e+06 ULP but the second patch does not improve anything.
I will post the first patch again so that that fix can be pushed to Newlib
but not the second one as it is only a MISRA coding standard fix. 
> Using RedHat newlib
> MPFR library: 4.1.0       
> MPFR header:  4.1.0 (based on 4.1.0)
> Checking function mylgammaf with MPFR_RNDN
> libm wrong by up to 7.50e+06 ulp(s) [7497618] for x=-0x1.3a7fcap+1
> mylgamma      gives -0x1p-24
> mpfr_mylgamma gives -0x1.e4cf24p-24
> Total: errors=509423944 (11.91%) errors2=11684280 maxerr=7.50e+06 ulp(s)
> 
> What result do you get for x=-0x1.3a7fcap+1?
> 
I get the same as you do.
> Best regards,
> Paul
> 
> > From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> > Date: Thu, 10 Feb 2022 17:12:40 +0100
> > 
> > The missing default case leads to large errors for |x| in range
> > [2.0, 3.0[.
> > In case of i=2, i.e., in the range between [2,3[, the computation
> > "r += logf(z)" is not performed and hence the result was wrong.
> > The added default case resolves this issue.
> > (Courtesy of Andreas Jung, ESA)
> > ---
> >  newlib/libm/math/er_lgamma.c  | 1 +
> >  newlib/libm/math/erf_lgamma.c | 1 +
> >  2 files changed, 2 insertions(+)
> > 
> > diff --git a/newlib/libm/math/er_lgamma.c b/newlib/libm/math/er_lgamma.c
> > index 5c88548fb..65727c6ab 100644
> > --- a/newlib/libm/math/er_lgamma.c
> > +++ b/newlib/libm/math/er_lgamma.c
> > @@ -302,6 +302,7 @@ static double zero=  0.00000000000000000000e+00;
> >            case 5: z *= (y+4.0);        /* FALLTHRU */
> >            case 4: z *= (y+3.0);        /* FALLTHRU */
> >            case 3: z *= (y+2.0);        /* FALLTHRU */
> > +      default:
> >                    r += __ieee754_log(z); break;
> >            }
> >      /* 8.0 <= x < 2**58 */
> > diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
> > index 84d02159b..e7311cacf 100644
> > --- a/newlib/libm/math/erf_lgamma.c
> > +++ b/newlib/libm/math/erf_lgamma.c
> > @@ -238,6 +238,7 @@ static float zero=  0.0000000000e+00;
> >            case 5: z *= (y+(float)4.0); /* FALLTHRU */
> >            case 4: z *= (y+(float)3.0); /* FALLTHRU */
> >            case 3: z *= (y+(float)2.0); /* FALLTHRU */
> > +      default:
> >                    r += __ieee754_logf(z); break;
> >            }
> >      /* 8.0 <= x < 2**58 */
> > -- 
> > 2.35.1

-- 
Andoni Arregi
Geschäftsführer
 
GTD GmbH
Ravensburger Str. 32a, 88677 Markdorf
T: +49 7544 96440 22 | M: +49 151 65620499 | F: +49 7544 96440 29
http://www.gtd-gmbh.de andoni.arregi@gtd-gmbh.de

Re: [PATCH 1/2] Improve lgammaf range for very small cases

[Andoni Arregi]

Hi Corinna,

On Thu, 2022-02-10 at 15:58 +0100, Corinna Vinschen wrote:
> Hi Andoni,
> 
> On Feb  9 17:28, Andoni Arregi wrote:
> > The original cut for small arguments at |x|<2**-70 (copied from the
> > double version) produces that when computing nadj we get a subnormal
> > number for t*x and thus, the division of pi/subnormal will be INF and
> > the logarithm of it too, which is wrong as a result for lgammaf in this
> > range.
> > The proposed new limit seems to be safe and has been tested to
> > produce accurate results.
> > (Courtesy of Andreas Jung, ESA)
> > ---
> >  newlib/libm/math/erf_lgamma.c | 2 +-
> >  1 file changed, 1 insertion(+), 1 deletion(-)
> > 
> > diff --git a/newlib/libm/math/erf_lgamma.c
> > b/newlib/libm/math/erf_lgamma.c
> 
> This is broken.  It should be on one line, so apparently your
> mailer broke it via autowrapping lines.
> 
Thanks for letting me know. Now trying it again without autowrapping.
> > index f88f63092..84d02159b 100644
> > --- a/newlib/libm/math/erf_lgamma.c
> > +++ b/newlib/libm/math/erf_lgamma.c
> > @@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
> >                *signgamp = -1;
> >            return one/(x-x);
> >         }
> > -       if(ix<0x1c800000) {     /* |x|<2**-70, return -log(|x|) */
> > +       if(ix<0x30800000) {     /* |x|<2**-30, return -log(|x|) */
> >            if(hx<0) {
> >                *signgamp = -1;
> >                return -__ieee754_logf(-x);
> > -- 
> > 2.35.1
> 
> But even if I fix this manually, your patch doesn't apply, neither with
> `git am', nor with `patch -p1'.  I have to admit that I don't understand
> why...
> 
I just checked `git am' myself before resubmitting. Hope it works now.
> Another question, would you mind to send patch series with a cover letter,
> please?
> 
Sure, sorry for that.
> 
> Thanks,
> Corinna

-- 
Andoni Arregi
Geschäftsführer
 
GTD GmbH
Ravensburger Str. 32a, 88677 Markdorf
T: +49 7544 96440 22 | M: +49 151 65620499 | F: +49 7544 96440 29
http://www.gtd-gmbh.de andoni.arregi@gtd-gmbh.de

Re: [PATCH 1/2] Improve lgammaf range for very small cases

[Andoni Arregi]

Hi Corina,

after my mistake from Tuesday and the checks by Paul Zimmermann
only the last patch I submitted for lgmammaf should be pushed as
Paul could confirm that the fix is correct. Forget about the
following patch series please:

[PATCH 0/2] Fix lgammaf and lgamma for small ranges and in [2, 3[ 

Regards

On Thu, 2022-02-10 at 15:58 +0100, Corinna Vinschen wrote:
> Hi Andoni,
> 
> On Feb  9 17:28, Andoni Arregi wrote:
> > The original cut for small arguments at |x|<2**-70 (copied from the
> > double version) produces that when computing nadj we get a subnormal
> > number for t*x and thus, the division of pi/subnormal will be INF and
> > the logarithm of it too, which is wrong as a result for lgammaf in this
> > range.
> > The proposed new limit seems to be safe and has been tested to
> > produce accurate results.
> > (Courtesy of Andreas Jung, ESA)
> > ---
> >  newlib/libm/math/erf_lgamma.c | 2 +-
> >  1 file changed, 1 insertion(+), 1 deletion(-)
> > 
> > diff --git a/newlib/libm/math/erf_lgamma.c
> > b/newlib/libm/math/erf_lgamma.c
> 
> This is broken.  It should be on one line, so apparently your
> mailer broke it via autowrapping lines.
> 
> > index f88f63092..84d02159b 100644
> > --- a/newlib/libm/math/erf_lgamma.c
> > +++ b/newlib/libm/math/erf_lgamma.c
> > @@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
> >                *signgamp = -1;
> >            return one/(x-x);
> >         }
> > -       if(ix<0x1c800000) {     /* |x|<2**-70, return -log(|x|) */
> > +       if(ix<0x30800000) {     /* |x|<2**-30, return -log(|x|) */
> >            if(hx<0) {
> >                *signgamp = -1;
> >                return -__ieee754_logf(-x);
> > -- 
> > 2.35.1
> 
> But even if I fix this manually, your patch doesn't apply, neither with
> `git am', nor with `patch -p1'.  I have to admit that I don't understand
> why...
> 
> Another question, would you mind to send patch series with a cover letter,
> please?
> 
> 
> Thanks,
> Corinna
> 

-- 
Andoni Arregi
Geschäftsführer
 
GTD GmbH
Ravensburger Str. 32a, 88677 Markdorf
T: +49 7544 96440 22 | M: +49 151 65620499 | F: +49 7544 96440 29
http://www.gtd-gmbh.de andoni.arregi@gtd-gmbh.de

[PATCH 2/2] Add a missing default case in lgamma

[Andoni Arregi]

The missing default case leads to large errors for |x| in range
[2.0, 3.0[.
In case of i=2, i.e., in the range between [2,3[, the computation
"r += logf(z)" is not performed and hence the result was wrong.
The added default case resolves this issue.
(Courtesy of Andreas Jung, ESA)
---
 newlib/libm/math/er_lgamma.c  | 1 +
 newlib/libm/math/erf_lgamma.c | 1 +
 2 files changed, 2 insertions(+)

diff --git a/newlib/libm/math/er_lgamma.c
b/newlib/libm/math/er_lgamma.c
index 5c88548fb..65727c6ab 100644
--- a/newlib/libm/math/er_lgamma.c
+++ b/newlib/libm/math/er_lgamma.c
@@ -302,6 +302,7 @@ static double zero=  0.00000000000000000000e+00;
 	   case 5: z *= (y+4.0);	/* FALLTHRU */
 	   case 4: z *= (y+3.0);	/* FALLTHRU */
 	   case 3: z *= (y+2.0);	/* FALLTHRU */
+      default:
 		   r += __ieee754_log(z); break;
 	   }
     /* 8.0 <= x < 2**58 */
diff --git a/newlib/libm/math/erf_lgamma.c
b/newlib/libm/math/erf_lgamma.c
index 84d02159b..e7311cacf 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -238,6 +238,7 @@ static float zero=  0.0000000000e+00;
 	   case 5: z *= (y+(float)4.0);	/* FALLTHRU */
 	   case 4: z *= (y+(float)3.0);	/* FALLTHRU */
 	   case 3: z *= (y+(float)2.0);	/* FALLTHRU */
+      default:
 		   r += __ieee754_logf(z); break;
 	   }
     /* 8.0 <= x < 2**58 */
-- 
2.35.1

Re: regression in Newlib 4.2.0 (lgammaf)

[Andoni Arregi]

Dear Paul,

here again, although we can provide a fix, we cannot understand well
enough why this didn't seem to be present in Newlib 4.1.0. The source
code regarding the incorrect limit at |x|<2**-70 has been there since
decades.

Best regards

On Sat, 2022-02-05 at 07:38 +0100, Paul Zimmermann wrote:
>        Hi,
> 
> a third regression in Newlib 4.2.0 for single precision:
> 
> Checking lgamma with newlib-4.2.0.20211231
> MPFR library: 4.1.0       
> MPFR header:  4.1.0 (based on 4.1.0)
> Checking function mylgammaf with MPFR_RNDN
> libm wrong by up to 8.93e+43 ulp(s) for x=-0x1p-70
> mylgamma      gives inf
> mpfr_mylgamma gives 0x1.842994p+5
> 
> The largest error for lgammaf in Newlib 4.1.0 was 7.50e6 ulps
> (for x=-0x1.3a7fcap+1).
> 
> Paul

-- 
Andoni Arregi
Geschäftsführer
 
GTD GmbH
Ravensburger Str. 32a, 88677 Markdorf
T: +49 7544 96440 22 | M: +49 151 65620499 | F: +49 7544 96440 29
http://www.gtd-gmbh.de andoni.arregi@gtd-gmbh.de

[PATCH] Improve lgammaf range for very small cases

[Andoni Arregi]

The original cut for small arguments at |x|<2**-70 (copied from the
double version) produces that when computing nadj we get a subnormal
number for t*x and thus, the division of pi/subnormal will be INF and
the logarithm of it too, which is wrong as a result for lgammaf in this
range.
The proposed new limit seems to be safe and has been tested to
produce accurate results.
(Courtesy of Andreas Jung, ESA)
---
 newlib/libm/math/erf_lgamma.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
index f88f63092..84d02159b 100644
--- a/newlib/libm/math/erf_lgamma.c
+++ b/newlib/libm/math/erf_lgamma.c
@@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
 	       *signgamp = -1;
 	   return one/(x-x);
 	}
-	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
+	if(ix<0x30800000) {	/* |x|<2**-30, return -log(|x|) */
 	   if(hx<0) {
 	       *signgamp = -1;
 	       return -__ieee754_logf(-x);
-- 
2.35.1

Re: [PATCH] Improve lgammaf range for very small cases

[Corinna Vinschen]

On Feb 11 12:16, Andoni Arregi wrote:
> The original cut for small arguments at |x|<2**-70 (copied from the
> double version) produces that when computing nadj we get a subnormal
> number for t*x and thus, the division of pi/subnormal will be INF and
> the logarithm of it too, which is wrong as a result for lgammaf in this
> range.
> The proposed new limit seems to be safe and has been tested to
> produce accurate results.
> (Courtesy of Andreas Jung, ESA)
> ---
>  newlib/libm/math/erf_lgamma.c | 2 +-
>  1 file changed, 1 insertion(+), 1 deletion(-)
> 
> diff --git a/newlib/libm/math/erf_lgamma.c b/newlib/libm/math/erf_lgamma.c
> index f88f63092..84d02159b 100644
> --- a/newlib/libm/math/erf_lgamma.c
> +++ b/newlib/libm/math/erf_lgamma.c
> @@ -168,7 +168,7 @@ static float zero=  0.0000000000e+00;
>  	       *signgamp = -1;
>  	   return one/(x-x);
>  	}
> -	if(ix<0x1c800000) {	/* |x|<2**-70, return -log(|x|) */
> +	if(ix<0x30800000) {	/* |x|<2**-30, return -log(|x|) */
>  	   if(hx<0) {
>  	       *signgamp = -1;
>  	       return -__ieee754_logf(-x);
> -- 
> 2.35.1
> 

Thanks, pushed.

However, there's still something wrong with your patch file.  In all
lines with spaces following tabs, your patch is missing exactly one
space.  Please check this locally.


Thanks,
Corinna