[Numpy-discussion] Different results from repeated calculation, part 2

[Xavier Gnata]

Alok Singhal wrote:
> On 14/08/08: 10:20, Keith Goodman wrote:
>   
>> A unit test is attached. It contains three tests:
>>
>> In test1, I construct matrices x and y and then repeatedly calculate z
>> = calc(x,y). The result z is the same every time. So this test passes.
>>
>> In test2, I construct matrices x and y each time before calculating z
>> = calc(x,y). Sometimes z is slightly different. But the x's test to be
>> equal and so do the y's. This test fails (on Debian Lenny, Core 2 Duo,
>> with libatlas3gf-sse2 but not with libatlas3gf-sse).
>>
>> test3 is the same as test2 but I calculate z like this: z =
>> calc(100*x,y) / (100 * 100). This test passes.
>>
>> I get:
>>
>> ======================================================================
>> FAIL: repeatability #2
>> ----------------------------------------------------------------------
>> Traceback (most recent call last):
>>   File "/home/[snip]/test/repeat_test.py", line 73, in test_repeat_2
>>     self.assert_(result, msg)
>> AssertionError: Max difference = 2.04946e-16
>>     
>
> Could this be because of how the calculations are done?  If the
> floating point numbers are stored in the cpu registers, in this case
> (intel core duo), they are 80-bit values, whereas 'double' precision
> is 64-bits.  Depending upon gcc's optimization settings, the amount of
> automatic variables, etc., it is entirely possible that the numbers
> are stored in registers only in some cases, and are in the RAM in
> other cases.  Thus, in your tests, sometimes some numbers get stored
> in the cpu registers, making the calculations with those values
> different from the case if they were not stored in the registers.
>
> See "The pitfalls of verifying floating-point computations" at
> http://portal.acm.org/citation.cfm?doid=1353445.1353446 (or if that
> needs subscription, you can download the PDF from
> http://arxiv.org/abs/cs/0701192).  The paper has a lot of examples of
> surprises like this.  Quote:
>
>   We shall discuss the following myths, among others:
>
>   ...
>
>   - "Arithmetic operations are deterministic; that is, if I do z=x+y in
>     two places in the same program and my program never touches x and y
>     in the meantime, then the results should be the same."
>   
>   - A variant: "If x < 1 tests true at one point, then x < 1 stays true
>     later if I never modify x."
>   
>   ...
>
> -Alok
>
>   

yep!
The code is buggy.
Please **nerver** use == to test if two floating point numbers are equal.
**nerver**.
As explained by Alok, floating point computations are *not* the same as 
computations over the field of reals numbers.

Intel 80bits registers versus 64bits representation, IEE754, sse or not 
sse or mmx : Fun with floating point arithmetic :)

The same code with and without SSE could provide you with "no equal" 
results.

Never use == on flaots?!? Well, you would have to use it in some corner 
cases but please remember that computers are only working on a small 
subset (finite) of natural numbers.

http://docs.python.org/tut/node16.html is a simple part of the story. 80 
bits register are the ugly part (but the fun one :))

abs(a-b)<epsilon is the correct way to test that a and b are "equal".

Xavier