Different results of the same caculation by loop code and by array code

r=1e9;
A=0;
for i=1:r
    A=A+floor(sqrt(r^2-i^2));
end
A

A = 7.8540e+17

---------------

i=[1:r]; 

Requested 1000000000x1 (7.5GB) array exceeds maximum array size preference (5.0GB). This might cause MATLAB to become unresponsive.

A=sum(floor(sqrt(r^2-i.^2)))

-------------------------

The same result by r=10^8

But the different results by r=10^9

Answer 1

[John D'Errico]

As is almost always the case, when we see this complaint, it is because you do not understand how arithmetic works in terms of double precision, and the capability of a double to correctly store numbers with wide dynamic ranges.

When r is 1e9, what is the SQUARE of r? How does that compare to the largest integer a double can store EXACTLY? What is flintmax for a double?

flintmax('double')

ans = 9.0072e+15

1e9^2

ans = 1.0000e+18

So by squaring r when it is that large, then performing arithmetic operations on it can cause the least significant digits of that large number to be incorrect. Looking more carefully at r^2, we see things like this

r = 1e9;
num2str(r^2,55)

ans = '1000000000000000000'

num2str(r^2 - 1^2,55)

ans = '1000000000000000000'

num2str(r^2 - 2^2,55)

ans = '1000000000000000000'

num2str(r^2 - 3^2,55)

ans = '1000000000000000000'

num2str(r^2 - 4^2,55)

ans = '1000000000000000000'

Do you see that in each case, subtracting a small number from r^2 did nothing? Now try the same thing using 1e8.

r = 1e8;
num2str(r^2,55)

ans = '10000000000000000'

num2str(r^2 - 1^2,55)

ans = '10000000000000000'

num2str(r^2 - 2^2,55)

ans = '9999999999999996'

num2str(r^2 - 3^2,55)

ans = '9999999999999992'

num2str(r^2 - 4^2,55)

ans = '9999999999999984'

In fact, even those results are failing sometimes down in the very last digit. This is because even 1e8^2 is greater than flintmax, but not by a lot.

In terms of getting a different result from a looped versus arrray code, you need to understand that performing the same computations in a subtly different sequence can yield a subtly different result. The classic example is:

0.1 + 0.2 - 0.3

ans = 5.5511e-17

-0.3 + 0.1 + 0.2

ans = 2.7756e-17

Neither computation yields zero, and the near zero result you see is different in the two mathematically equivalent cases.

And so we see tiny errors down in the least significant bits. But when you compare a looped versus an array computation for the "same" thing, you need to understand that MATLAB will control what sequence the vectorized computation is done in, and that is done because you want to see your computations done quickly. In turn, that means farming out those computations to multiple cores in your CPU. So the sum in the end need not happen in eactly the same order. Again, two computations done in a different sequence on numbers can and often will result in a different result.

The end result is a tradeoff between speed, which surely you demand, and the least significant bits of the result.

[dpb]

Does MATLAB also perhaps do something like Kahan summation inside sum as opposed to explicit linear looping code?

[Walter Roberson]

No, MATLAB does not do Kahan summation.

vectorized summation, for large enough arrays, is done by splitting the vector into pieces, using optimized C++ code to form partial sums on each of the pieces, and then add the partial sums together.