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

[Daniel J Sebald]

On 12/15/2014 09:53 AM, Jordi Gutiérrez Hermoso wrote:
> On Sun, 2014-12-14 at 23:05 -0800, s.jo wrote:
> > I found some errors in my Cygwin system:
>
> Hardly an error,
>
>      http://wiki.octave.org/FAQ#Why_is_this_floating_point_computation_wrong.3F
>
> if you are writing code that depends on checking exact zeros instead
> of almost zeros, your code is buggy. Yes, your code may have a bug.
> You cannot depend on exact equality except in very narrow cases.
>
> > Can other users check this example to see if all cygwin release has same
> > problem?
>
> For your case, I cannot reproduce what you consider a problem:
>
>      octave:1>  a = randn
>      a = -1.6148
>      octave:2>  b = randn
>      b = -0.73621
>      octave:3>  a^2 + b^2
>      ans =  3.1497
>      octave:4>  a^2 + b^2 == a*a + b*b
>      ans =  1
>      octave:5>  [a,b; -b a]*[a;b]
>      ans =
>
>         3.14966
>         0.00000
>
>      octave:6>  [a,b; -b a]*[a;b] == [a^2 + b^2; 0]
>      ans =
>
>         1
>         1
>
> The fact that you're seeing errors of the order of 1e-18, which is
> smaller than eps(1), makes me think that your BLAS compilation is
> probably storing some floats in registers, giving you a different
> precision.

This is probably correct.  The expression a*a + b*b is mathematically 
the same as the 1,1 element of the matrix multiplication, but that 
doesn't necessarily mean the machine code is compiled in exactly the 
same way.  For example, in the matrix multiplication case it might be 
that there is some cumulative sum (after all, generally the matrix 
operation works on N pairs).

I would guess though that it is possible to make the underlying code the 
same in the parser/interpreter somehow with a little work.  The cost 
would likely be a small loss of efficiency given we're looking at a very 
fundamental feature.  Is it worth it?  Not sure, but it's probably a 
more logical thing to attempt than trying to declare equality when 
floats are next to one another.

Dan