Date: Fri, 8 Mar 2019 08:24:32 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@freebsd.org Subject: Re: cexpl() (was: Re: Update ENTERI() macro) Message-ID: <20190308162432.GA31392@troutmask.apl.washington.edu> In-Reply-To: <20190308225807.D6410@besplex.bde.org> References: <20190306055201.GA40298@troutmask.apl.washington.edu> <20190306225811.P2731@besplex.bde.org> <20190306184829.GA44023@troutmask.apl.washington.edu> <20190307061214.R4911@besplex.bde.org> <20190306214233.GA23159@troutmask.apl.washington.edu> <20190307144315.N932@besplex.bde.org> <20190307044447.GA16298@troutmask.apl.washington.edu> <20190307163958.V1293@besplex.bde.org> <20190307195436.GA20856@troutmask.apl.washington.edu> <20190308225807.D6410@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Sat, Mar 09, 2019 at 12:51:07AM +1100, Bruce Evans wrote:
> On Thu, 7 Mar 2019, Steve Kargl wrote:
>
> > On Thu, Mar 07, 2019 at 05:36:22PM +1100, Bruce Evans wrote:
> >> On Wed, 6 Mar 2019, Steve Kargl wrote:
> >>
> >>> On Thu, Mar 07, 2019 at 03:09:06PM +1100, Bruce Evans wrote:
> >* ...
> >>> /*
> >>> * exp_ovfl = 11356.5234062941439497, 0x400c, 0xb17217f7d1cf79ac
> >>> * cexp_ovfl = 22755.3287906024445633, 0x400d, 0xb1c6a8573de9768c
> >>> */
>
> I checked that these are correct. Any rounding should be down for
> exp_ovfl and up for cexp_ovfl.
>
> >>> if ((hx == 0x400c && lx > 0xb17217f7d1cf79acULL) ||
> >>> (hx == 0x400d && lx < 0xb1c6a8573de9768cULL)) {
> >>
> >> Don't duplicate the hex constants in the comment. The limits can be fuzzy
> >> so don't need so many decimal digits or nonzero nonzero bits or LD80C()
> >> to round them (s_cexp.c uses only 1 decimal digit and doesn't test the
> >> lower 32 bits in the mantissa; here we have more bits than we need in lx).
> >
> > Being fuzzy was the movitation for the new macro, I was only using
> > the upper 32 bits of lx in the original cexpl() I posted.
>
> IIRC, that was done by extracting into a 32-bit lx. This could be written
> as (lx >> 32) > 0xb17217f7 (or a fuzzier value), and it might be better
> to do that anyway.
This is what I do in my WIP ccosh. There are 3 thresholds
of ~22.9, ~11356, and ~22755. I'm casting the result to
uint32_t, so something like (uint32_t)(lx >> 32) > 0xb17000000.
Do I need the cast and/or do you prefer it?
> >> Maybe not use bit-fiddling at all for this (in other precisions too).
> >> The more important exp() functions check thresholds in FP. Here x can
> >> be NaN so FP comparisons with it don't work right. That seems to be
> >> the only reason to check in bits here.
> >>
> >>> /*
> >>> * x is between exp_ovfl and cexp_ovfl, so we must scale to
> >>> * avoid overflow in exp(x).
> >>> */
> >>> RETURNI(__ldexp_cexpl(z, 1));
> >>
> >> Did your tests cover this case, and does it use the fixed version of the
> >> old __ldexp_cexpl() in k_expl.h?
> >
> > No. The test I reported was restricted to -11350 < x < 11350.
> > Spot checked x outside that range. I had the conditional wrong
> > for x < -exp_ovfl. Also, the comment is actaully a little optimistic.
> > One cannot avoid overflow. For example, x = 11357 is just above
> > exp_ovfl, and exp(x) = inf. __ldexp_cexpl() does scaling and exploits
> > |cos(y)| <= 1 and |sin(y)| <= 1 to try to extend the range of cexpl().
>
> I understand this better now. __ldexp_cexpl() doesn't work for all large
> x, so callers must only use it for a range of values above cexp_ovl.
> The correctness of this is unobvious. It depends on sin(y) and cos(y)
> never underflowing to 0. Otherwise, this method would give Inf * 0 = NaN
> for finite x and y, and the scaling method would give finite * 0 * 2**k = 0.
cos(y) can never be 0 and sin(y) is only 0 at y = 0 (see sinpi, cospi
implementations). y=0 is a special case and handled separately. We
can have Inf*subnormal when y is near 0.
--
Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20190308162432.GA31392>