[Numpy-discussion] Different results from repeated calculation, part 2
Fri Aug 15 20:02:52 CDT 2008
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."
The code is buggy.
Please **nerver** use == to test if two floating point numbers are equal.
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"
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".
More information about the Numpy-discussion