Re: What to do about tgammal?

From: Steve Kargl <sgk_at_troutmask.apl.washington.edu>
Date: Tue, 14 Dec 2021 23:14:22 UTC
On Tue, Dec 14, 2021 at 03:47:17PM -0700, Warner Losh wrote:
> On Tue, Dec 14, 2021 at 3:23 PM Mark Murray <markm@freebsd.org> wrote:
> 
> > On 14 Dec 2021, at 21:51, Steve Kargl <sgk@troutmask.apl.washington.edu>
> > wrote:
> > > Interval          max ULP          x at Max ULP
> > > [6,1755.1]        0.873414 at 1.480588145237629047468e+03
> > > [1.0662,6]        0.861508 at 1.999467927053585410537e+00
> > > [1.01e-17,1.0661] 0.938041 at 1.023286481537296307856e+00
> > > [-1.9999,-1.0001] 3.157770 at -1.246957268051453610329e+00
> > > [-2.9999,-2.0001] 2.987659 at -2.220949465449893090070e+00
> > >
> > > Note, 1.01e-17 can be reduced to soemthing like 1.01e-19 or
> >
> > Extra diffs most welcome!
> >
> 
> Those results have allayed my fears. Is 1e-17 sufficient to commit these
> changes, or do
> we need to work to get it down to 1e-19 for it to be acceptable?
> 

Well, egg-on-the-face.  I updated the sloppy limit to isolate
x = 0 from 0x1p-56 to use 0x1p-116 and forgot to update my test
scripts.  The code in tgammal in the interval [-iota:iota] is

	if (x > iota)
		RETURNI(smaller_gam(x));

	if (x > -iota) {
		if (x != 0)
			u.a = 1 - tiny;	/* raise inexact */
		RETURNI(1 / x);
	}

so tgammal(x) -> 1/x as x -> 0.


-- 
Steve