Re: What to do about tgammal?
- In reply to: Warner Losh : "Re: What to do about tgammal?"
- Go to: [ bottom of page ] [ top of archives ] [ this month ]
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