Re: ***[Possible UCE]*** Re: More efficient "dbl_raise"

[Daniel J S. <dan...@ie...>]

Hans-Bernhard Bröker wrote:
> Daniel J Sebald wrote:
> 
> 
>>The issue in dbl_raise is a multiplying loop for raising a value to an integer 
>>power.  I'm guessing the library function pow() does such a thing more 
>>efficiently in addition to handling the pathological NaN.  
> 
> 
> Maybe so, but at the same time, there's a solid probability that it is 
> just too inaccurate to work.
> 
> Tics generation is extremely sensitive to rounding error, mainly in 
> extreme cases: very large or small arguments, arguments already slightly 
> distrubed by earlier rounding errors.
> 
> A while ago we had a long-standing bug where tic generation was broken 
> just by turning on -O2 in GCC, but only on some platforms.  Suffice it 
> to say that this is a kind of problem I'd rather not have to face again.
> 
> 
>>In some way it also removes the question that Hans raised in the code.  I'm sure 
>>you can't recall that far back, Hans, but a question I have is just exactly was 
>>this approach supposed to achieve?  
> 
> 
> Basically, it's that pow() isn't required to treat integer exponents 
> specially, so it can easily destroy more precision than dbl_raise() in 
> its existing form does.

Right.  But the premise, I guess, is that for the application in question we're 
interested in raising an integer to an integer power, which results in an 
integer.  So, we are free to round the pow() result to an integer and can feel 
fairly certain the result is the same as dbl_raise().

Your statement may be true for relatively small numbers.  To confirm that for 
larger numbers, one would have to look at the pow() code.  I mean, the way 
gnuplot is programmed, what if someone enters numbers on the order of 1e30? 
What's the precision of 10*10*...*10*10 as a floating point number?  It's beyond 
the largest value the mantissa can hold so we know some precision is likely lost 
since the underlying math is binary floating point.

Also, there may even be something about the existing code that is questionable, 
but probably the log10() function behaves nicely in this regard.  This bit:

floor(log10(arg))

We aren't guaranteed that will come out to be exactly the exponent desired.  Say 
the argument is 1e15 for which we'd expect the resulting value to be 15.  But, 
if by chance because of the way the particular log10() function is implemented 
in a library, it comes out to 14.9999992, the floor() would produce the 
incorrect result.  We can't round in this case.

This kind of precision issue probably works best if the number system is binary 
float.  Then we could use frexp() and ldexp() and get the exponent directly and 
do highly accurates tests and checks.  But I don't see any creative way of 
changing the problem into that realm.

Maybe the best we could hope for is some kind of internal consistency between 
log10() and pow(), by which I mean it would be nice that if

   E = floor(log10(arg))

then

   pow(10,E) <= arg < pow(10,E+1)

is always true.  For what it is worth, this extra test:

exponent = floor(log10(arg));
if (pow(10,exponent+1) <= arg)
     exponent++;

does nothing to change the PostScript results.

Dan