Date: Wed, 6 Mar 2019 13:42:33 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@freebsd.org Subject: Re: Update ENTERI() macro Message-ID: <20190306214233.GA23159@troutmask.apl.washington.edu> In-Reply-To: <20190307061214.R4911@besplex.bde.org> References: <20190227074811.GA75972@troutmask.apl.washington.edu> <20190227201214.V1823@besplex.bde.org> <20190227161906.GA77785@troutmask.apl.washington.edu> <20190228060920.R4413@besplex.bde.org> <20190304212159.GA12587@troutmask.apl.washington.edu> <20190305153243.Y1349@besplex.bde.org> <20190306055201.GA40298@troutmask.apl.washington.edu> <20190306225811.P2731@besplex.bde.org> <20190306184829.GA44023@troutmask.apl.washington.edu> <20190307061214.R4911@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20190306214233.GA23159>