[Matrix-SIG] possible bug in NumPy

[Tim Peters]

[Oliver Gathmann]
> ...
> I think I found what looks like a bug in NumPy (I am not up
> to date with the very latest release, but I don't remember
> this coming up here before):

I haven't yet made time to get NumPy at all, but believe I can shed some
painful light regardless <wink>.

> Python 1.5 (#3, Feb 12 1998, 12:14:31)  [GCC 2.7.2.3] on linux2
> Copyright 1991-1995 Stichting Mathematisch Centrum, Amsterdam
> >>> from Numeric import *
> >>> a = array([[1,1],[2,2]])
> >>> b = array([[1.1,1.1],[2.1,2.1]])
> >>> c = array([[0.1,0.1],[0.1,0.1]])
> >>> abs(a-b)
> array([[0.1,0.1],
>        [0.1,0.1]])
> >>> less_equal(abs(a-b),c)
> array([[0, 0],
>        [0, 0]])
> >>> less_equal(c,abs(a-b))
> array([[1, 1],
>        [1, 1]])
>
> Shouldn't that be 'all true' in both cases since abs(a-b) and c are
> the same thing?

*If* they were the same thing, yes.  But forget arrays for a moment:

Python 1.5.1 (#0, Apr 13 1998, 20:22:04) [MSC 32 bit (Intel)] on win32
Copyright 1991-1995 Stichting Mathematisch Centrum, Amsterdam
>>> x = 1.1 - 1
>>> y = .1
>>> x <= y, y <= x
(0, 1)

I.e., the same thing happens for scalars.

>>> x == y
0
>>> x - y
8.32667268469e-017

I.e., floating-point roundoff caught you in the last bit -- neither .1 nor
1.1 are exactly representable as doubles, so y isn't actually one tenth, and
while "1.1 - 1" suffers no roundoff error, the double version of "1.1"
wasn't actually one and one tenth to begin with either (and loses
approximately the last 3 bits of the approximation to .1 in order to "make
room for" the leading 1).

Clearer?

I don't know what you should *do* about it in NumPy, but someone else here
will -- this kind of problem has existed since the first floating-point
computation ever done <wink>.

what-you-see-isn't-always-what-you-get-ly y'rs  - tim