Update ENTERI() macro
Bruce Evans
brde at optusnet.com.au
Wed Mar 6 19:30:49 UTC 2019
On Wed, 6 Mar 2019, Steve Kargl wrote:
> On Wed, Mar 06, 2019 at 11:56:23PM +1100, Bruce Evans wrote:
>> On Tue, 5 Mar 2019, Steve Kargl wrote:
>
>>> a similar k_cexpl.c. Yes, I added the 'c' in the name to avoid
>>> confusion in ld80/. In particular, I have no idea how he found
>>> his scaling value 'k'. Any insights?
>>
>> bde already wrote __ldexp_cexpl() in ld*/k_expl.h, and kargl committed
>> it in r260066. Does it not work? :-).
>>
>> Well, it hasn't been tested, and it indeed cannot work since it spells
>> cosl(y) as cos(y).
>
> Taking long breaks from pecking at libm issues seems to be
> conducive to memory loss. I'll go review k_expl.h. I simply
> remember it as having a kernel for expl().
I now see that you implemented 2 more versions of __ldexpl_cexpl() by
cloning the old double precision version. Apparently the includes
are to unpolluted for the compiler to see the multiple versions :-).
Using the version in k_expl.h almost forces inlining of expl()'s kernel
and its large tables, just like for hyperbolic functions. This wastes
a lot of space, especially for duplicating the tables. It is only a
small optimization for time. It is done for the hyperbolic functions
to get this optimization, and for __ldexpl_cexpl() just for convenience.
>> ...
>> I have rewritten the double and float versions in work related to
>> fixing the accuracy of the double and float versions of hyperbolic
>> functions. I fixed these before writing the long double hyperbolic
>> ..
>> XX Index: k_exp.c
>
> Any chance you'll get around to committing your WIP? Yes,
> I know you have a few thousand kernel patches in your
> queue above libm patches. :-)
Probably not soon.
>
> BTW, for the non-exceptional cases with 1M random z values
> where z=x+Iy and -11350 < x,y < 11350 and I only consider
> results that are normal, I find my cexpl() yields
>
> % ./testl -u -X 11350
> Max ULP Re: 1.980535
> z = (1.22918109220546510585e+03,1.03853909865862542237e+04)
> libm = (5.06780736805320166327e+533,-4.39418989082799451477e+533)
> mpfr = (5.06780736805320166270e+533,-4.39418989082799451456e+533)
> Max ULP Im: 2.022155
> z = (4.83490728165160559637e+03,1.07778990242305355345e+04)
> libm = (-3.66535128319537945953e+2099,4.67021177841072936494e+2099)
> mpfr = (-3.66535128319537945915e+2099,4.67021177841072936441e+2099)
Check a few denormals, Infs and NaNs.
I can only easily test for coherence with double precision. That is,
evaluate in both precisions and check that the long double results
rounded down are within a few ulps.
> For comparison, cexp() with -705 < x,y < 705 yields
>
> % ./testd -u -X 705
> Max ULP Re: 2.215132
> z = (1.49377521822925502e+02,1.79997882645095970e+01)
> libm = (4.93720465697268180e+64,-5.61754869856313932e+64)
> mpfr = (4.93720465697268310e+64,-5.61754869856313987e+64)
> Max ULP Im: 2.182779
> z = (2.50664219501672335e+02,-4.81327697040560906e+02)
> libm = (-5.73256778461974670e+108,4.48612518733245315e+108)
> mpfr = (-5.73256778461974754e+108,4.48612518733245429e+108)
>
> Certainly, not exhaustive but encouraging.
We expect long double precision to be slightly more accurate, at
least using the expl kernel, since the kernel accuracy is about
0.51 ulps while the exp accuracy is only about 0.8 ulps. The
above shows 0.2 ulps better instead of 0.3. Still more than 2
ulps total from combining several errors of 0.5-1.0 ulps.
Bruce
More information about the freebsd-numerics
mailing list