Re: odd behaviour with 2x2 matrix floating-point multiplication

[s.jo]

After Jordi's comment, I updated my cygwin system from 1.7.32
to 1.7.33-2 (ver  2014-11-13).

Now I have a perfect computation on zero entry, but a*a+b*b term still shows 
mismatch with a probability of 0.24 or so: If I try 10000 pairs of a and b,
about 2300 cases shows mismatch at the a*a+b*b.
I checked that the mismatch occurs at only last bit of hex form, say
+/-eps(a*a+b*b),
as Dan pointed out.

I understand that N pairs term floating-point operation of matrix could be
different
from the exact math expression.
But, in this case, only two terms are involved in machine code that BLAS may
generate. 
I am not sure of how BLAS routine can shuffle the operation order or do
anything
to show different floating-point results.
I suspect that BLAS routines scale the numbers unnecessarily, even trivial
cases.

Is this a problem related to a release version of BLAS and/or cygwin?

-- jo

ps.

The follow script is to check the probability of a*a+b*b mismatch:

%% 
num_trial=10000;
count_err=[0;0];
for ii=1:num_trial
        a = randn;
        b = randn;
        C = [a b; -b a]*[a; b];
        D = [a*a+b*b;-b*a+a*b];
        err=(C ~= D);
        count_err = count_err + err;
        if count_err(1)==100
                num2hex(C(1))
                num2hex(D(1))
                num2hex(D(1)-eps(D(1)))
                num2hex(D(1)+eps(D(1)))
                break
        end
end
ii
count_err
ii/count_err

%% result at 100th mismatch 
ans = 40181cd9a6c0611d
ans = 40181cd9a6c0611e
ans = 40181cd9a6c0611d
ans = 40181cd9a6c0611f
a =    -1.96601051483154
b =     1.47070618400453
ans = 40181cd9a6c0611e
ii =                  415
count_err =

                   100
                     0
ans = 

     0.240963855421687
                     0




--
View this message in context: 
http://octave.1599824.n4.nabble.com/odd-behaviour-with-2x2-matrix-floating-point-multiplication-tp4667799p4667812.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.