small cleanup patch for e_pow.c

[Bruce Evans]

On Sun, 10 May 2015, Steve Kargl wrote:

> On Mon, May 11, 2015 at 06:36:51AM +1000, Bruce Evans wrote:
>> On Sun, 10 May 2015, Steve Kargl wrote:
>>>
>>> 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.
>
> I used p = 24 in my (1+-2^(-p)).  It seems that you're using p-1=23.

It should be 1+2**-(p-1) and 1-2**-p, but the 1-epsilon side dominates;
also 2**128 is sloppy with gradual underflow; in the scale factor for
the 1-epsilon side, 128 must be replaced by something like 150.  So
30.7 is correct if we don't bother to treat the sides differently.

There is a related problem with extra precision.  Expressions like
huge*huge and tiny*tiny are always wrong in the presence of extra
precision/exponent range, since they don't overflow or underflow
but just give a wrong result if extra exponent range actually works.
But extra precision/exponent is often broken:
- C11 broke it and also broke efficiency by requiring extra precision/
   exponent range to be destroyed on function return.  clang doesn't
   implement this breakage, but gcc-4.8 does (with -std=c99 but not
   with -std=gnu99)
- at least old versions of clang miscompile expressions like huge*huge
   unless huge is declared with the volatile hack.  The C11 breakage
   makes returning the value INFINITY correct and then the only error
   in evaluating this at compile time is not raising the overflow
   exception.
- at least old versons of both clang and gcc miscompile expressions
   like tiny*tiny, except...

Functions like exp() are especially affected by this bug.  cexp() is
most noticeably broken by it, since it uses expressions like
exp(x) * cos(y), written otherwise carefully by scaling exp(x) so
as to avoid spurious overflow, so some non-overflowing cases return
garbage related to huge*huge in extra precision.

We mostly don't worry about this bug, except I fixed it for cexp*().
My test don't see it, since they assign function results to a variable
without the extra precision, so as to destroy it like the C11 breakage
requires.  The bug should be fixed someday by replacing code like
'return huge*huge;' with 'return Inf_with_OE();'.  Cases where the
expression is actually x*x where x is known to be >= huge are harder
to fix.

Bruce