small cleanup patch for e_pow.c

[Steve Kargl]

On Mon, May 11, 2015 at 06:36:51AM +1000, Bruce Evans wrote:
> 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.
> 

I used p = 24 in my (1+-2^(-p)).  It seems that you're using p-1=23.

(1-2**(-p))**(2**N) = 2**(emin-p)

Using log() and log(1+z) = z for |z| << 1, the above equation gives

N = p + log((p-emin) * log(2)) / log(2)

where the equality is of course approximation.

p = 24 --> N = 30.690
p = 23 --> N = 29.681

-- 
Steve