Re: More efficient "dbl_raise"

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

Ethan A Merritt wrote:
> On Saturday 07 April 2007 21:28, you wrote:
> 
>>Ethan A Merritt wrote:
>>
>>>On Saturday 07 April 2007 15:58, Daniel J Sebald wrote:
>>>
>>>
>>>>Move the patch into CVS or dump it?  (Discussing things and not resolving them,
>>>>even for as trivial a patch as this, isn't productive.)
>>>
>>>
>>>Was there a real problem that this was supposed to solve?
>>>If not, drop it.
>>
>>Recap:  The patch mimics dbl_raise() behavior, but uses an existing library 
>>function pow() that is probably more efficient
> 
> 
> I was not paying close attention, but didn't HBB show that your
> proposed "more efficient" alternative function could be off by
> 15 orders of magnitude  in the worst case?

Well, I'm not sure.  This might be an apples and oranges kind of thing.  Here is 
the script HBB suggested:

     powdiff(base,power) = (base**power / exp(power*log(base))) - 1.0
     set samples 301 ; plot [0:60] powdiff(10,x)

Note that this sampling contains both "integers" and "floats".  What we are 
ultimately interested in our use of pow(,) or dbl_raise() is raising an integer 
value (speaking math-wise, not computer-wise) to an integer power (again, 
math-wise), but that is getting off course.  Let's pause to look at the behavior 
of gnuplot on some simple examples, say 10^10:

gnuplot> print 10**10
1410065408
gnuplot> print 10**10.0
10000000000.0
gnuplot> print 10.0**10
10000000000.0
gnuplot> print exp(10*log(10))
10000000000.0

pow(base,power) = base**power

gnuplot> print pow(10,10)
1410065408
gnuplot> print pow(10,10.0)
10000000000.0

OK, so we see that gnuplot treats integers (computer-wise) with integer math 
which is much more limiting in size compared to floating point math.

But this is not how the innards of the pow(), floor(), etc. functions work:

#include <math.h>
   double floor( double arg );

The documentation of these C functions might say "largest integer not greater 
than arg", but that means the math vernacular, not computer vernacular.  That is 
why I say apples and oranges.

Given gnuplot's behavior, can we come up with a comparison at the command line 
the reflects the behavior of C routines?  I've followed the parsing to 
internal.c and see that ** operator does in fact use the C library pow(,) 
function.  So we've got that.  The expression exp(power*log(base)) loses 
accuracy, but that isn't really related to the use of pow() I've applied in the 
range application.  So let's not look at that any further.

Try this:

powint(base,power) = power==0 ? 1.0 : base * powint(base,power - 1)
powdiff(base,power) = (base**power / powint(base,power)) - 1.0
set samples 61 ; plot [0:60] powdiff(10,x)

You may see something different than what I'm seeing, but what I see is that 
near x=25 there appears to be some discrepancy.  So let's pick a few points there:

gnuplot> print 10.0**23 - powint(10.0,23)
0.0
gnuplot> print 10.0**24 - powint(10.0,24)
0.0
gnuplot> print 10.0**25 - powint(10.0,25)
2147483648.0
gnuplot> print 10.0**26 - powint(10.0,26)
17179869184.0
gnuplot> print 10.0**27 - powint(10.0,27)
137438953472.0
gnuplot> print 10.0**28 - powint(10.0,28)
0.0
gnuplot> print 10.0**29 - powint(10.0,29)
0.0

So, something does become flaky here.  Let's see if we can figure out where 
exactly the discrepancy is:

gnuplot> print 10.0**23 - 1e23
0.0
gnuplot> print 10.0**24 - 1e24
0.0
gnuplot> print 10.0**25 - 1e25
0.0
gnuplot> print 10.0**26 - 1e26
0.0
gnuplot> print 10.0**27 - 1e27
0.0
gnuplot> print 10.0**28 - 1e28
0.0
gnuplot> print powint(10.0,23.0) - 1e23
0.0
gnuplot> print powint(10.0,24.0) - 1e24
0.0
gnuplot> print powint(10.0,25.0) - 1e25
-2147483648.0
gnuplot> print powint(10.0,26.0) - 1e26
-17179869184.0
gnuplot> print powint(10.0,27.0) - 1e27
-137438953472.0
gnuplot> print powint(10.0,28.0) - 1e28
0.0
gnuplot>

Now that' odd.  It would seem that the power function is correct and multiplying 
isn't.  Let's go one step further to confirm this:

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
1e+25
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

Very strange.  Not concrete proof, but my conjecture about not knowing which is 
more accurate, outright multiplying via dbl_raise() vs. pow(), does seem a 
pertinent question.

I was long winded, Ethan, but the answer to your question is "I don't think so."

Dan