(0.3)^3 == (0.3)*(0.3)*(0.3) returns false in matlab?

Question

I am trying to understand roundoff error for basic arithmetic operations in MATLAB and I came across the following curious example.

(0.3)^3 == (0.3)*(0.3)*(0.3)

ans = 0

I'd like to know exactly how the left-hand side is computed. MATLAB documentation suggests that for integer powers an 'exponentiation by squaring' algorithm is used.

"Matrix power. X^p is X to the power p, if p is a scalar. If p is an integer, the power is computed by repeated squaring."

So I assumed (0.3)^3 and (0.3)*(0.3)^2 would return the same value. But this is not the case. How do I explain the difference in roundoff error?

Answer 1

I don't know anything about MATLAB, but I tried it in Ruby:

irb> 0.3 ** 3
  => 0.026999999999999996
irb> 0.3 * 0.3 * 0.3
  => 0.027

According to the Ruby source code, the exponentiation operator casts the right-hand operand to a float if the left-hand operand is a float, and then calls the standard C function pow(). The float variant of the pow() function must implement a more complex algorithm for handling non-integer exponents, which would use operations that result in roundoff error. Maybe MATLAB works similarly.

Answer 2

Interestingly, scalar ^ seems to be implemented using pow while matrix ^ is implemented using square-and-multiply. To wit:

octave:13> format hex
octave:14> 0.3^3
ans = 3f9ba5e353f7ced8
octave:15> 0.3*0.3*0.3
ans = 3f9ba5e353f7ced9
octave:20> [0.3 0;0 0.3]^3
ans =

  3f9ba5e353f7ced9  0000000000000000
  0000000000000000  3f9ba5e353f7ced9

octave:21> [0.3 0;0 0.3] * [0.3 0;0 0.3] * [0.3 0;0 0.3]
ans =

  3f9ba5e353f7ced9  0000000000000000
  0000000000000000  3f9ba5e353f7ced9

This is confirmed by running octave under gdb and setting a breakpoint in pow.

The same is likely true in matlab, but I can't really verify.

Answer 3

Thanks to @Dougal I found this:

#include <stdio.h>
int main() {
  double x = 0.3;
  printf("%.40f\n", (x*x*x));

  long double y = 0.3;
  printf("%.40f\n", (double)(y*y*y));
}

which gives:

0.0269999999999999996946886682280819513835
0.0269999999999999962252417162744677625597

The case is strange because the computation with more digits gives a worst result. This is due to the fact that anyway the initial number 0.3 is approximated with few digits and hence we start with a relatively "large" error. In this particular case what happens is that the computation with few digits gives another "large" error but with opposite sign... hence compensating the initial one. Instead the computation with more digits gives a second smaller error but the first one remains.

Answer 4

Here's a little test program that follows what the system pow() from Source/Intel/xmm_power.c, in Apple's Libm-2026, does in this case:

#include <stdio.h>
int main() {
    // basically lines 1130-1157 of xmm_power.c, modified a bit to remove
    // irrelevant things

    double x = .3;
    int i = 3;

    //calculate ix = f**i
    long double ix = 1.0, lx = (long double) x;

    //calculate x**i by doing lots of multiplication
    int mask = 1;

    //for each of the bits set in i, multiply ix by x**(2**bit_position)
    while(i != 0)
    {
        if( i & mask )
        {
            ix *= lx;
            i -= mask;
        }
        mask += mask;
        lx *= lx; // In double this might overflow spuriously, but not in long double
    }

    printf("%.40f\n", (double) ix);
}

This prints out 0.0269999999999999962252417162744677625597, which agrees with the results I get for .3 ^ 3 in Matlab and .3 ** 3 in Python (and we know the latter just calls this code). By contrast, .3 * .3 * .3 for me gets 0.0269999999999999996946886682280819513835, which is the same thing that you get if you just ask to print out 0.027 to that many decimal places and so is presumably the closest double.

So there's the algorithm. We could track out exactly what value is set at each step, but it's not too surprising that it would round to a very slightly smaller number given a different algorithm for doing it.