Re: What to do about tgammal?

From: Mark Murray <markm_at_FreeBSD.org>
Date: Sat, 18 Dec 2021 10:41:14 UTC

> On 18 Dec 2021, at 03:52, Steve Kargl <sgk@troutmask.apl.washington.edu> wrote:
> 
> On Tue, Dec 14, 2021 at 10:21:55PM +0000, Mark Murray 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!
>> 
> 
> Hi Mark,
> 
> Don't know if you noticed, but I borroewed a few cpu cycles
> from grimoire.

Didn't notice a thing! :-)

>  My WIP is already better than the imprecise.c
> kludge from theraven@.  I need to work out the details of
> computing logl(x) in extra precision or see if I can leverage
> what Bruce did a few years ago.  Anywho, current results:
> 
> Interval tested for tgammal: [128,1755]
> count: 1000000
>  xm =  1.71195767195767195767195767195767183e+03L
> libm =  7.79438030237108165017007606176403036e+4790L
> mpfr =  7.79438030237108165017007606175285456e+4790L
> ULP = 14869.19517
> 
> Interval tested for tgammal: [16,128]
> count: 1000000
>  xm =  1.27687183687183687183687183687183690e+02L
> libm =  6.61421998891483212224382625339007663e+212L
> mpfr =  6.61421998891483212224382625338960267e+212L
> ULP = 731.00958
> 
> Interval tested for tgammal: [10,16]
> count: 1000000
>  xm =  1.54261654261654261654261654261654251e+01L
> libm =  2.74203137295418912508367515208072654e+11L
> mpfr =  2.74203137295418912508367515208073861e+11L
> ULP = 45.61161
> 
> Interval tested for tgammal: [1.2446e-60,10]
> count: 1000000
>  xm =  6.26200626138006138006138006138006065e-02L
> libm =  1.54507103764516989381203274093299079e+01L
> mpfr =  1.54507103764516989381203274093299091e+01L
> ULP = 0.76751

Hmm. I think my understanding of ULP is missing something?

I thought that ULP could not be greater than the mantissa size
in bits?

I.e., I thought it represents average rounding error (compared with
"perfect rounding"), not truncation error, as the above very large
ULPs suggest.

M
--
Mark R V Murray