Update ENTERI() macro

Steve Kargl sgk at troutmask.apl.washington.edu
Wed Mar 6 18:48:35 UTC 2019

On Wed, Mar 06, 2019 at 11:56:23PM +1100, Bruce Evans wrote:
> On Tue, 5 Mar 2019, Steve Kargl wrote:
> > ...
> > I think I have the s_cexpl.c file fixed to use LDBL_EXTRACT_WORDS
> > instead of the new macro I introduced.  I however cannot figure
> > out what David Das did to arrive at k_exp.c, so I cannot write
> Davvid A S.

Whoops, again.  Seems I conflated real name and login name.

> > 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().

> It is written in a more portable way than the double and float versions
> using '1' instead of bit fiddling to start the construction of 2**k.
> This makes it identical in ld80 and ld128 (since the exponent range
> is the same and the accessor macros hide enough details of the bit
> fiddling for the exponent).  The duplication is noted in a comment,
> but the comment also says that using scalbnl() ond ld128 is probably
> best after all (it is slow, but multiplication is also slow).
> 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
> functions using identical methods.  You committed the latter, but
> the double and float versions still use the old innaccurate slower
> methods.
> The rewrite improves the comments, and it is the improved comments
> that the reference to k_exp.c in k_expl.h is supposed to be for.
> 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. :-)

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)

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.


More information about the freebsd-numerics mailing list