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

[Bruce Evans]

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.

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

>>> 	} else {
>>> 		/*
>>> 		 * Cases covered here:
>>> 		 *  -  x < exp_ovfl and exp(x) won't overflow (common case)
>>> 		 *  -  x > cexp_ovfl, so exp(x) * s overflows for all s > 0
>>> 		 *  -  x = +-Inf (generated by exp())
>>> 		 *  -  x = NaN (spurious inexact exception from y)
>>> 		 */
>>> 		exp_x = expl(x);
>>> 		sincosl(y, &s, &c);
>>> 		RETURNI(CMPLXL(exp_x * c, exp_x * s));
>>> 	}
>>> }
>>
>> I think the last clause can be simplfied and even optimized by using the
>> kernel in all cases:

I was wrong.

>> - the kernel produces exp_x in the form kexp_x * 2**k where kexp_x is near
>>    1.  Just multiply kexp_x by c and s, and later multiply by 2**k.  This
>>    cost an extra multiplication by 2**k (for the real and imaginary parts
>>    separately), but avoids other work so may be faster.
>> - hopefully kexp_x is actually >= 1 (and < 2).  Otherwise we have to
>>    worry about spurious underflow
>> - the above method seems to avoid spurious underflow automatically, and
>>    doesn't even mention this (multiplying by exp_x may underflow, but not
>>    spuriously since c and s are <= 1).
>
> You might be right about using the kernel.  I'll leave that
> to others on freebsd-numerics@; although you and I seem to be
> the only subscribers.

Forget about that.  __ldexp_cexp*() is not really a kernel.  It is used
for ccosh*() and csinh*() too.  Those can only use the kernel for large
positive x.  cexp*() may as well do the same.  Also, the kernel currently
doesn't support very large positive x, so the second overflow threshold
is needed in callers.

Bruce