computing expontential as repeated multiplication for small integer exponent



>>>>> "Robert" == Robert Dodier <robert.dodier at gmail.com> writes:

    Robert> On 7/30/07, Raymond Toy <raymond.toy at ericsson.com> wrote:
    >> Your example of 1.2345678901234567d0^3 just illustrates that expt and
    >> multiplication are different. As best as I can tell, (expt x 3) from
    >> cmucl (and sbcl?) is more accurate than repeated multiplication.

    Robert> At this point, I think it is more important to get various Lisp
    Robert> implementations to agree, rather than getting maximally-accurate
    Robert> results for some implementations. I say this after having crunched

I prefer maximally-accurate results if possible.  I'm more than
capable of getting somewhat accurate results, so I want (someone else)
to give me maximally-accurate results. :-)

    Robert> more than a few numbers in various contexts.

    Robert> That said, I would be equally happy with an implementation of
    Robert> EXPTB which didn't treat small integer exponents as repeated
    Robert> multiplication. All I want is for different Lisps to get the same result.

Why is this so important?  What if sin(1.23456789101112d0) returns
different results from different lisps?  Are you going to add a
floating-point sin implementation to maxima just so they all agree?
And what about the Lisp reader?  There are lots of floating-point
coefficients in the Lisp source.  Are you going to force all lisps to
use the same reader to produce the same numbers?

    >> I think rtest8 #69 shouldn't expect exactly the same result.  We should
    >> add a fuzz factor, like in some of the other tests for approximately equal.

    Robert> Well, this fuzz factor business is problematic. There is already
    Robert> a fuzz factor; it is 8 * (64 bit float epsilon). Test #69 comes out

Yes, that the fuzz factor varies is a problem.

    Robert> with a greater error because it is the result of an iterative algorithm
    Robert> (namely find_root); a small discrepancy in EXPTB causes find_root
    Robert> to take a different number of steps for different Lisps.

    Robert> In general it is difficult to gauge an acceptable discrepancy.
    Robert> For one of the test cases for lsquares which I was recently working
    Robert> on, different Lisps agree to only 3 decimal places.

But doesn't this mean that either your problem is ill-conditioned,
or your algorithm is unstable?  Do you know the analytical solution?
If not how to you know which one is "correct"?

    Robert> When floating point operations aren't guaranteed to come
    Robert> out the same, we can't find a suitable fuzz factor without

But we can't guarantee that FP operations will come out the same
unless maxima implements them all.  I don't think we want to do that.
The best we can do is assume IEEE arithmetic and assume the 4 basic
ops plus sqrt are correct.  

Don't get me wrong.  I understand your desire to have all Lisps
produce the same result.  I just don't think making things less
accurate is the way to go.  (I have asserted that sbcl has an accurate
expt.  I haven't proven it, though.)

Ray