MATLAB: Cumsum gives different results for double and single

cumsumcumsum;double;single;weirddoublesingle

Hey,
Does anyone have an explanation for this?
a = single(rand(1000,1000));
b = cumsum(cumsum(a));
c = single(cumsum(cumsum(double(a))));
sum(b(:)~=c(:))
ans =
      864947
Thanks!
Yuval

Best Answer

single() and double() are handled by hardware instructions using IEEE 754 standards on floating point arithmetic.
Suppose you add
1/100000 + 1/99999 + 1/99998 + .... 1/3 + 1/2 + 1
and you use a particular number of bits in the mantissa. As you add terms, the cumulative sum becomes large relative to the values of the terms being added, and some terms may effectively give no contribution -- as the "focus" of the floating point number shifts "left" as the sum gets larger, some terms may be less than 1/2 the contribution of the least significant bit of the sum, and so will not contribute anything to the sum. They might as well not be there.
Now run the same sum again, with the same original approximations, but with more bits in the mantissa. 53 bits to accumulate into takes a lot longer to overflow than 24 bits, so those bits that would formerly be lost from contributing start to contribute -- and they add up! For example the terms between 1/10000 and 1/1001 add up to about 2.3 . The terms between 2^(-53) and 2^(-24) add up to about 20.1, and that is a contribution that can be meaningfully made at 53 bits that would be lost at 24 bits
You are mentally picturing an addition scheme in which no bits are lost as you add those single precision numbers, sort of like adding integers and using as many bits as you need. But that is not how adding floating point numbers works: any part of the next value that does not fit within the accumulating value is just thrown away, but the contributions of those parts you are throwing away can turn out to be important and would be more retained with a wider register.
Related Question