* [PATCH 0/2] Fix lgammaf and lgamma for small ranges and in [2, 3[
2022-02-10 14:58 ` Corinna Vinschen
@ 2022-02-10 16:07 ` Andoni Arregi
2022-02-10 16:10 ` [PATCH 1/2] Improve lgammaf range for very small cases Andoni Arregi
` (3 subsequent siblings)
4 siblings, 0 replies; 15+ messages in thread
From: Andoni Arregi @ 2022-02-10 16:07 UTC (permalink / raw)
To: newlib
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
^ permalink raw reply [flat|nested] 15+ messages in thread
* [PATCH 1/2] Improve lgammaf range for very small cases
2022-02-10 14:58 ` Corinna Vinschen
2022-02-10 16:07 ` [PATCH 0/2] Fix lgammaf and lgamma for small ranges and in [2, 3[ Andoni Arregi
@ 2022-02-10 16:10 ` Andoni Arregi
2022-02-10 16:21 ` Paul Zimmermann
2022-02-10 16:12 ` [PATCH 2/2] Add a missing default case in lgamma Andoni Arregi
` (2 subsequent siblings)
4 siblings, 1 reply; 15+ messages in thread
From: Andoni Arregi @ 2022-02-10 16:10 UTC (permalink / raw)
To: newlib
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
^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: [PATCH 1/2] Improve lgammaf range for very small cases
2022-02-10 16:10 ` [PATCH 1/2] Improve lgammaf range for very small cases Andoni Arregi
@ 2022-02-10 16:21 ` Paul Zimmermann
0 siblings, 0 replies; 15+ messages in thread
From: Paul Zimmermann @ 2022-02-10 16:21 UTC (permalink / raw)
To: andoni.arregui; +Cc: newlib
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
>
>
>
^ permalink raw reply [flat|nested] 15+ messages in thread
* [PATCH 2/2] Add a missing default case in lgamma
2022-02-10 14:58 ` Corinna Vinschen
2022-02-10 16:07 ` [PATCH 0/2] Fix lgammaf and lgamma for small ranges and in [2, 3[ Andoni Arregi
2022-02-10 16:10 ` [PATCH 1/2] Improve lgammaf range for very small cases Andoni Arregi
@ 2022-02-10 16:12 ` Andoni Arregi
2022-02-10 17:29 ` Paul Zimmermann
2022-02-10 16:16 ` [PATCH 1/2] Improve lgammaf range for very small cases Andoni Arregi
2022-02-11 11:19 ` Andoni Arregi
4 siblings, 1 reply; 15+ messages in thread
From: Andoni Arregi @ 2022-02-10 16:12 UTC (permalink / raw)
To: newlib
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
^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: [PATCH 2/2] Add a missing default case in lgamma
2022-02-10 16:12 ` [PATCH 2/2] Add a missing default case in lgamma Andoni Arregi
@ 2022-02-10 17:29 ` Paul Zimmermann
2022-02-11 11:14 ` Andoni Arregi
0 siblings, 1 reply; 15+ messages in thread
From: Paul Zimmermann @ 2022-02-10 17:29 UTC (permalink / raw)
To: andoni.arregui; +Cc: newlib
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
^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: [PATCH 2/2] Add a missing default case in lgamma
2022-02-10 17:29 ` Paul Zimmermann
@ 2022-02-11 11:14 ` Andoni Arregi
0 siblings, 0 replies; 15+ messages in thread
From: Andoni Arregi @ 2022-02-11 11:14 UTC (permalink / raw)
To: Paul Zimmermann; +Cc: newlib
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
^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: [PATCH 1/2] Improve lgammaf range for very small cases
2022-02-10 14:58 ` Corinna Vinschen
` (2 preceding siblings ...)
2022-02-10 16:12 ` [PATCH 2/2] Add a missing default case in lgamma Andoni Arregi
@ 2022-02-10 16:16 ` Andoni Arregi
2022-02-11 11:19 ` Andoni Arregi
4 siblings, 0 replies; 15+ messages in thread
From: Andoni Arregi @ 2022-02-10 16:16 UTC (permalink / raw)
To: newlib
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
^ permalink raw reply [flat|nested] 15+ messages in thread
* Re: [PATCH 1/2] Improve lgammaf range for very small cases
2022-02-10 14:58 ` Corinna Vinschen
` (3 preceding siblings ...)
2022-02-10 16:16 ` [PATCH 1/2] Improve lgammaf range for very small cases Andoni Arregi
@ 2022-02-11 11:19 ` Andoni Arregi
4 siblings, 0 replies; 15+ messages in thread
From: Andoni Arregi @ 2022-02-11 11:19 UTC (permalink / raw)
To: newlib
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
^ permalink raw reply [flat|nested] 15+ messages in thread