Re: More efficient "dbl_raise"

[Ethan A M. <merritt@u.washington.edu>]

On Sunday 08 April 2007 01:52, Daniel J Sebald wrote:
> Daniel J Sebald wrote:
> 
> > gnuplot> print 
> > 10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0*10.0 
> > - 1e25
> > -2147483648.0
> 
> And this seems to be accurate:
> 
> gnuplot> print 
> 1e25/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10/10 
> - 1.0
> 0.0
> 
> Ummmm, so why is the division accurate while the multiplication above appears 
> not to be?  

That is exactly what you should expect.  Sequential division reduces the inaccuracy
of the original number by a factor of 10 at each step.  Simple subtraction exposes
the size of the original discrepancy.

Note that 2147483648.0, although a large number in absolute terms, corresponds
to an error in the 15th decimal place of the decimal representation of 10^25.
More to the point, that corresponds to about 50 bits of precision.
Since double-precision IEEE format only provides 52 bits of precision in the
mantissa, this is within one bit of the expected rounding error.
(And that one bit may be my mistake; I'm doing this on my fingers, as it were).

Note also that if you care about this level of precision then you also must
pay attention to the IEEE rounding mode.

> Is there something going on with the translation from exponential  
> notation and the machine IEEE float representation somewhere around 10^25?

Yes. That's where you hit the limit of the IEEE representation.
To do better than that you'd have to switch to some other floating point
representation.

There are infinite-precision math libraries available, but gnuplot isn't
using one.

-- 
Ethan A Merritt