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

[Andoni Arregi <andoni.arregui@gtd-gmbh.de>]

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