Update ENTERI() macro

[Steve Kargl]

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:
> >> On Wed, 6 Mar 2019, Steve Kargl wrote:
> >> ...
> >> But does it preserve NaN bits as much as possible, and produce the same
> >> NaN bits as double precision after conversion to double precision?  The
> >> report spells NaN and Inf using C99's fuzzy bad names.  Even the C99
> >> specification spells NaN and Inf correctly in Appendix F.
> >>
> >> The complex hyperbolic functions have better comments, so are almost
> >> an easier place to start.
> >
> > I suppose it preserves the NaN bits as well as cexpf and cexp as
> > it uses the same algorthim.
> 
> > ...
> > long double complex
> > cexpl (long double complex z)
> > {
> > 	long double c, exp_x, s, x, y;
> > 	uint64_t lx, ly;
> > 	uint16_t hx, hy, ix, iy;
> >
> > 	ENTERI();
> >
> > 	x = creall(z);
> > 	y = cimagl(z);
> >
> > 	EXTRACT_LDBL80_WORDS(hx, lx, x);
> > 	EXTRACT_LDBL80_WORDS(hy, ly, y);
> 
> I think you recently changed to this from your new macro.

Yes, the new macro is gone.

> >
> > 	ix = hx&0x7fff;
> > 	iy = hy&0x7fff;
> 
> Some style bugs.  s_cexp.c doesn't use the fdlibm'ish style of sometimes
> no spaces around '&'.  It puts the ix and iy calculations immediately
> after the extractions.  They shouldn't be separated since they are
> also extractions.  c_cexp.c extracts in a different order, and doesn't
> use iy.  You reordered it a bit to use sincos().  I don't like the
> reordering much.  But old s_cexp.c has its own style bug of a separating
> the extraction for y but not even following its usual style of a blank
> line before comments near the extraction of x.

I may unorder the re-ordering to the old order.  If I understand
one of the comments in the GCC bugzilla report, the optimization of 
converting

  a = cos(x);
  b = sin(x):

to 

  sincos(x, &b, &a).

Actually, does "a = cos(x); b = sin(x)" goes to cexp(cmplx(0,x)) in
the middle-end of the compiler.  The backend then can recognize
cexp(cmplx(0,x)) and converts this to sincos(x, &b, &a).  I can't 
verify this for FreeBSD, because libm isn't C99 compliant.  GCC
requires a C99 compliant to activate the optimization. 

I fixed the extraction for hx, and hy.

> > 	/* cexp(0 + I y) = cos(y) + I sin(y) */
> > 	if ((ix | lx) == 0) {
> > 		sincosl(y, &s, &c);
> > 		RETURNI(CMPLXL(c, s));
> > 	}
> >
> > 	/* cexp(x + I 0) = exp(x) + I 0 */
> > 	if ((iy | ly) == 0)
> > 		RETURNI(CMPLXL(expl(x), y));
> >
> > 	if (iy == 0x7fff) {
> 
> s_cexp.c actually uses a clobbered hy here.  Perhaps fix this too (compilers
> mostly ignore program order for initializing variables like iy and optimize
> to an expression involving the unclobbered hy here if that is optimal.  The
> extra variables are just for readability or unreadability).
> 
> > 		if (ix < 0x7fff ||
> > 		    (ix == 0x7fff && (lx & 0x7fffffffffffffffULL) != 0)) {
> > 			/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
> > 			RETURNI(CMPLXL(y - y, y - y));
> 
> Looks OK.
> 
> I think ix and iy are only used here, so it is clearer to not have
> variables for them -- just use (hx & 0x7fff), etc.
> 
> y - y and similar expressions usually give the right NaN (y = Inf -> dNaN
> and y = NaN -> d(y)).
> 
> IIRC, the only differences for ld128 is that then normalization bit is
> hidden so masking out the top bit in the mantissa is not needed, but the
> mantissa is in 2 64-bit words.
> 
> > 		} else if (hx & 0x8000) {
> 
> This is clearest as it is with an explicit mask.
> 
> > 	/*
> > 	 *  exp_ovfl = 11356.5234062941439497, 0x400c, 0xb17217f7d1cf79ac
> > 	 * cexp_ovfl = 22755.3287906024445633, 0x400d, 0xb1c6a8573de9768c
> > 	 */
> > 	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.  

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

 ./testl -a 11357 0.1
libm: Re(cexpf) =  inf
mpfr: Re(cexpf) =  1.90658369228334884925e+4932
libm: Im(cexpf) =  1.91296449568717353572e+4931
mpfr: Im(cexpf) =  1.91296449568717353558e+4931

troutmask:kargl[430] ./testl -a 11357 1.57
libm: Re(cexpf) =  1.52588659766018377153e+4929
mpfr: Re(cexpf) =  1.52588659766018377147e+4929
libm: Im(cexpf) =  inf
mpfr: Im(cexpf) =  1.91615588587371861485e+4932


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

-- 
Steve