cexpl() (was: Re: Update ENTERI() macro)

[Bruce Evans]

On Fri, 8 Mar 2019, Steve Kargl wrote:

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

I think you can even round to a power of 2 and not look at any mantissa
bits to classify this, as for the 64 in coshl().

22.9 here seems wrong.  It is approximately the threshold where cosh()
decays to exp().  The threshold for coshl() delay is much higher.
coshl() uses 64, which is high enough for both ld80 and ld128.  More
like 32 would work for ld80.

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

coshl() also checks thresholds in only in FP if classifications using the
exponent are too fuzzy, so as to avoid the slow step of extracting the
mantissa bits and the unportable steps of using them.  I think
classification for ld80 has to look at the mantissa to classify UNnormals,
but ld128 doesn't have UnNormals so might not need the mantissa bits.

Look at coshl() for details.  It is simpler.
>>
>> 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.

cos(y) and sin(y) would be zero if the exact value underflows on
assignment to long double, or the inexact value is miscalculated as
0.  It is difficult to prove that they don't underflow to 0 or to a
denormal when y is not near 0.  That they are not miscalculated so as
to give spurious underflow follows if we can prove that denormals don't
even nearly occur and that our methods have errors less than 1 ulp in
relative terms (only values at least twice as large as the largest
denormal occur, and a 1-ulp relative error for these values can't reach
a denormal value).

Inf*denormal is not very special.  It is +-Inf.  The upper overflow threshold
is chosen so than exp(x) * +-smallest_denormal = +-Inf for x above this
threshold.

It is easier to optimize this for cexp().  cexp() starts with extracting
all bits of y, and then all bits of x unless y = 0, but in most cases the
args are finite, nonzero and not very large, and everything except zero
can be classified by looking at only the exponent.  Zero can be classified
using an FP comparison.

Something like:

 	/* Do this first for gcc's sin+cos pessimization? */
 	if (x == 0)
 		return (cpackl(cosl(y), sinl(y));

 	if (y == 0)
 		return (cpackl(expl(x), y));

 	/*
 	 * Maybe do the following earlier and check (hx & 0x7fff) == 0,
 	 * etc., before the FP classification of 0.
 	 */
 	hx = GET_LDBL_EXPSIGN(x);
 	hy = GET_LDBL_EXPSIGN(y);

 	if (hx & (0x7fff) < MUMBLE && (hy & 0x7fff) != 0x7fff)
 		goto final_clause_now;	/* simple formula */

 	/* Need more bits for a fuller classification (rarely reached). */
 	lx = GET_LDBL_MANTISSA(x);
 	ly = GET_LDBL_MANTISSA(y);

This is supposed to handle UnNormals without explicit checks by letting
other layers and hardware handle them however it wants.  E.g., if x or
y is Pseudo-Zero, then the comparision with 0 either succeeds or fails
depending on whether Pseudo-Zero is treated as 0 or NaN by the hardware
comparison operator, and our treatment is consistent for the early returns.
We only have to worry about falling through to later code which is
inconsistent because it depends on Pseudo-Zeros not getting that far.

Bruce