Update ENTERI() macro

[Steve Kargl]

On Thu, Mar 07, 2019 at 06:30:42AM +1100, Bruce Evans wrote:
> 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.

The version in k_expl.h has 2 bugs.  You note the first (cos instead
of cosl).  The second is

In file included from /data/kargl/trunk/math/libm/msun/ld80/s_cexpl.c:43:
/data/kargl/trunk/math/libm/msun/ld80/k_expl.h:288:22: error: magnitude of
      floating-point constant too large for type 'double'; maximum is
      1.7976931348623157E+308 [-Werror,-Wliteral-range]
        exp_x = (lo + hi) * 0x1p16382;
                            ^
1 error generated.
*** Error code 1


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

Exceptional cases give

% ./testl -e
cexpl(1, inf) = (nan,nan) expecting (nan,nan)
cexpl(1,-inf) = (nan,nan) expecting (nan,nan)
cexpl(nan,-inf) = (nan,nan) expecting (nan,nan)
cexpl(nan, inf) = (nan,nan) expecting (nan,nan)
cexpl(inf, inf) = (inf,nan) expecting (inf,nan)
cexpl(inf,-inf) = (inf,nan) expecting (inf,nan)
cexpl(inf,nan) = (inf,nan) expecting (inf,nan)
cexpl(-inf,nan) = (0,0) expecting (0,0)
cexpl(-inf,-inf) = (0,0) expecting (0,0)
cexpl(-inf, inf) = (0,0) expecting (0,0)

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

I use GNU MPFR and compute exp(x)*cos(y) and exp(x)*sin(y)
with a precision of 4*LDBL_MANT_DIG.  Supposedly, each 
function is correctly rounded to this precision.  I take
the answers as-if it is exact when computing ULP.

> 
> > 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

-- 
Steve