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

       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 <[email protected]>
> 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
> 
> 
>