Differences calculating the sum of squares in single between sum() and a loop

Update: Please Visit the Following Simplified Version of This Question http://www.mathworks.com/matlabcentral/answers/85211-different-results-between-2013a-32bit-and-64bit-in-single-precision

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

Hi, I'm calculating the sum of squares in C and MATLAB in single precision and noted some differences. I wondered if anyone can fill me in on the behaviour below. Suspect floating point associative problem, but still seems a little odd...

Numbers in my array a are 'reasonably' well scaled, from a column in a single precision matrix generated by MATLAB. Differences being seen here are common and leading to substantial differences between same algorithm in C and MATLAB.

Answer displayed from code below:

ans =

   4.0864482
   4.0864477
   4.0864482
   4.0864482

Code to reproduce:

clc
format long
a = single([-0.112033270,0.153425455,-0.751111090,-0.129932076,0.206365064,-0.871111095,0.000000000,...
            0.000000000,0.000000000,0.000662919,-0.000005372,0.004444445,0.002651675,-0.000085950,...
            0.017777778,0.005966269,-0.000435120,0.039999999,0.010606701,-0.001375194,0.071111113,...
            0.016572969,-0.003357407,0.111111112,0.023865076,-0.006961920,0.159999996,0.032483019,...
            -0.012897816,0.217777774,0.042426802,-0.022003105,0.284444451,0.053696420,-0.035244718,...
            0.360000014,0.066291876,-0.053718518,0.444444448,0.080213174,-0.078649282,0.537777781,...
            0.095460303,-0.111390717,0.639999986,0.112033270,-0.153425455,0.751111090,0.129932076,...
            -0.206365064,0.871111095]);

%Method one (matlab code)
r1 = sum(a.^2);

%Method two (as done in C)
r2=single(0);
for i = 1:length(a)
    r2 = r2 + a(i)^2;
end

%Method 3 (variation, correct answer (why?))
r3=single(0);
for i = 1:length(a)
    x = a(i)^2;
    r3 = r3 + x;
end

%Method 4 (why does double(0) change things - single takes precedence?)
r4 = double(0);
for i = 1:length(a)
    r4 = r4 + a(i)^2;
end

display([r1;r2;r3;r4])

[Roger Stafford]

Single precision floating point numbers have only 24 bits in their significands (mantissas), so you should expect errors of the order of 1 part in 2^24 or 2^23, or about seven decimal place accuracy. That's roughly what I see in your results. Why does that surprise you?

[Iain]

Looking at his code, it is not clear why the results of the operation would differ.

His results indicate that:

 r2 = r2 + a(i)^2;

and

    x = a(i)^2;
    r3 = r3 + x;

are different despite the fact that the operations are identical.

[Roger Stafford]

Those results could be different if in the first case the computer does the squaring and adding all in double precision before rounding back to single precision, while in the second case such a rounding might occur both after the squaring operation and after the addition.

These single precision results differ only in their` least bit. It takes only one rounding difference among the fifty-one steps to produce that.

[Jonathan Currie]

OK I have simplified the problem to two elements. New simplified question posed: http://www.mathworks.com/matlabcentral/answers/85211-different-results-between-2013a-32bit-and-64bit-in-single-precision

Answer 1

[the cyclist]

What version of MATLAB are you using? I get

ans =

     4.0864482
     4.0864482
     4.0864482
     4.0864482

from your code, using R2013a on a OS X 10.8.4.

[Jonathan Currie]

Hi - yep just checked on another PC, I get the correct answer on Matlab 2013a 64bit, but the incorrect answer on 2013a 32bit (Win64)... Can anyone else repeat this result?