small cleanup patch for e_pow.c

[Bruce Evans]

On Sun, 10 May 2015, Steve Kargl wrote:

> On Sun, May 10, 2015 at 08:16:14PM +1000, Bruce Evans wrote:
>> ...
>> So 1 is the only numbers near 1 that doesn't give overflow.
>
> Thanks for the explanation!  That help dislodge a mental block.
>
> To find the magic numbers, it seems I need to consider
>
> (1-2**(-p))**(2**N) = 2**(emin-p)          for underflow
> (1+2**(-p))**(2**N) = (1-2**(-p))*2**emax  for overflow
>
> With p = [24, 53, 64, 113], emin = [-125, -1021, -16381, -16381],
> emax = [128, 1024, 16384, 16384], and the use of log(1+z) = z for
> |z| << 1, I find
>
> underflow:  N = [30.7, 62.5, 77.5, 126.5]
> overflow:   N = [30.5, 62.5, 77.5, 126.5]

I plugged some numbers into pari to get sloppy estimates.  E.g.,
(1 + 2^-23.0)^(2^N) / 2^128; then bump up N until the result is > 1.
I get this at N = 29.5, so 30.7 seems too high.

> PS: Yes, I'm slowly poking away at powl().  I have the 19 special cases
> working for ld80.  Now, comes the hard part.

Excellent.

The only thing I would change immediately in the fdlibm implementation is
its manually inlined exp().  Use the exp*() kernel instead.  Hopefully
the kernel used for hyperbolic functions is also suitable for pow().
This kernel has only been committed for long double precision.  Look at
the old BSD libm for use of other exp() kernels and their use by pow().
It has 2 kernels: exp__D(), which is used by pow(), erf() and [t]gamma(),
and exp__E(), which is used by cosh() and expm1().  We still have
exp__D() (renamed with more underscores) in bsdsrc for use in tgamma().
A general kernel would take an extra-precision arg and return an extra-
precision result, but that would be inefficient and these kernels only
support an extra-precision arg.  I don't see how that is useful for
expm1(), and my kernel for exp*() only gives extra precision for the
result, as needed for hyperbolic functions.

I didn't notice before that old BSD libm used an exp() kernel for
cosh().  Ng (or other Sun hackers?) apparently dumbed down his
implementation of cosh() a bit for the version in fdlibm, so FreeBSD's
cosh() is probably still less accurate than the old (1985) version
:-(.  This is fixed in my version, but only the long double version
of the fixes has been committed.  The accuracy of cosh() was presumably
better than 1 ulp in 1985; it is 2 ulps (maybe 1.5) in FreeBSD, and
about 0.9 ulps in my version (the kernel for expl() is more accurate,
so coshl( has an accuracy of closer to 0.5 ulps).  OTOH, pow() may
have been improved over the 1985 version in fdlibm.  It manually inlines
kernels for both exp() and log() (named exp__D() and log__L(), and
these kernels seem to be accurate enough -- I haven't found any errors
larger than 1 ulp in pow().  The manual inlining makes the code harder
to read (about 50% larger, with many more special constants for the
kernels).

Bruce