Re: More efficient "dbl_raise"

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

Hans-Bernhard Bröker wrote:
> Daniel J Sebald wrote:
> 
>> Hans-Bernhard Bröker wrote:
>>
>>> Daniel J Sebald wrote:
> 
> 
>>>> 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.  
> 
> 
> No. We're interested in raising a double to an integer power, which 
> results in a double.  For simple cases, the input and output doubles 
> will have integer values, but they're not integers by design.

This equation

     double power = dbl_raise(10.0, floor(log10(arg)));

from a purely mathematical standpoint, is raising an integer, 10, to an integer, 
floor(), power.  Generally, that's a rational number.  Now, looking more closely 
at dbl_raise(x,y), we have

     int i = abs(y);

An integer raised to a positive integer power is an integer.  (I forgot to 
clarify it is a positive integer power last time, maybe that's where the 
confusion is.)

We are free to round the value and remove any finite math effects.


 >  If we
 > only had to worry about cases where the power is small enough for the
 > result to be integer, a table of all 10 powers-of-ten from 10^0 to 10^9
 > would suffice.  But that's not the case.
[snip]
 > As good as it can be.  It's still likely to be better than that of the
 > typical naive implementation of pow(x,y):

If the ultimate result can't be represented with a binary float, that's a 
different matter.  But I'm saying that isn't unique to pow() because neither can

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

be contained in binary float.  It reaches a point along the way in the 
multiplication loop that the ability to represent the large integer is gone and 
each multiplication results in poorer and poorer resolution relative to the 
increasing product.  Who knows?  Perhaps the pow() function has better numerical 
behavior and ends up being more accurate in some circumstances.


> For the fun of it, be sure to make some plots of this function:
> 
>     powdiff(base,power) = (base**power / exp(power*log(base))) - 1.0
> 
> E.g.
> 
>     set samples 301 ; plot [0:60] powdiff(10,x)
> 
> to see just how badly wrong this can go.
> 
>> floor(log10(arg))

Interesting plot, and that is the issue at hand.


>>
>> We aren't guaranteed that will come out to be exactly the exponent 
>> desired.  
> 
> 
> No, we're not.  Which is why this result is used only as a guide, not as 
> the single piece of information, by quantize_tics.  There's a reason 
> that there are cases of that switch outside the expected range of [2:20].
> 
>> 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.  
> 
> 
> No, the best we can do is avoid having to rely on such assumptions.  The 
>  code in quantize_tics() does that.

You're looking at this in one level broader scope than I am.  I'm just focused 
on getting xnorm to be accurate given whatever the value of arg is, i.e., the 
correct value of the exponent E rather than possibly being off by one... which 
is a little subjective in itself meaning that E being one less than the 
mathematically correct value won't result in overly bad number of tics anyway.

Dan