Python issues with rounding

Question

I have a rounding issue in Python 2.7 leading to unexpected output. I'm trying to get combinations of p1 and p2 that sum up to 0.6 or less in total.

from itertools import product
P = []
p1 = [0.0,0.2,0.4,0.6]
p2 = [0.0,0.2,0.4,0.6] 
for p11,p22 in product(p1,p2):
    if p11+p22 <= max(p1):    
        P.append((p11,p22))

However, when I run this, it does not include all values for which p11+p22 = 0.6:

[(0.0, 0.0),
 (0.0, 0.2),
 (0.0, 0.4),
 (0.0, 0.6),
 (0.2, 0.0),
 (0.2, 0.2),
 (0.4, 0.0),
 (0.6, 0.0)]

It does work correctly when I set p11+p22 <= max(p1)+0.01. For different p1 and p2 the problem may or may not occur. I find this behavior extremely strange, leading to very unreliable results.

It is probably related to floating precision issues. In my opinion this behaviour shouldn't exist in Python, since R and Matlab do not have this behaviour either. Are there any simple ways around this?

Answer 1

What is happening?

Computers have an internal representation of numbers. In most cases, those representations have a fixed number of bits. This leads to only a fixed amount of numbers being representable. For example, you might know that languages Like C have a maximum value for integers.

Similar, you can't store the exact representation of some floating point numbers. As the computer uses base two, there are some numbers in base 10 which have a short, finite representation but the binary one is long. For more details, see IEEE 754.

How can it be "fixed"?

There is nothing to be fixed here as everything is working like it was specified. But you have to know about these types of problems. When you are aware of the fact that there is a problem, then there are two strategies to get around it.

Either use epsilons (-> don't compare with exact numbers, but check if the number is within a very small interval around the number. The length of this interval is often called "epsilon") or use arbitrary precision representations (see fractions. The second only works when you can influence how the number is put into the program, e.g.

from itertools import product
from fractions import Fraction
P = []
p1 = [Fraction(0.0), Fraction(2, 10), Fraction(4, 10), Fraction(6, 10)]
p2 = [Fraction(0.0), Fraction(2, 10), Fraction(4, 10), Fraction(6, 10)]
for p11, p22 in product(p1, p2):
    if p11+p22 <= max(p1):
        P.append((p11, p22))

See also

• Is floating point math broken?
• What Every Computer Scientist Should Know About Floating-Point Arithmetic

Answer 2

Because the limitations on fixed bit width floating point, you need to use an arbitrary precision floating point package or explicitly compare with +/- an epsilon amount.

Python includes decimal (for 'arithmetic that works in the same way as the arithmetic that people learn at school'):

from itertools import product
import decimal

P = []
p1 = map(decimal.Decimal, ['0.0','0.2','0.4','0.6'])
p2 = map(decimal.Decimal, ['0.0','0.2','0.4','0.6'])
for p11,p22 in product(p1,p2):
    if p11+p22 <= max(p1):    
        P.append((p11,p22))

Answer 3

If you want to compare with matlab or R, or have performance issues , here is a numpy approach, with np.isclose() as workaround.

p1 = [0.0,0.2,0.4,0.6]
p2 = [0.0,0.2,0.4,0.6] 
sums=np.add.outer(p1,p2)
P1,P2=np.meshgrid(p1,p2)
compare = (sums<0.6) | np.isclose(sums,0.6)
print(np.dstack((P1,P2))[compare])

which gives :

[[ 0.   0. ]
 [ 0.2  0. ]
 [ 0.4  0. ]
 [ 0.6  0. ]
 [ 0.   0.2]
 [ 0.2  0.2]
 [ 0.4  0.2]
 [ 0.   0.4]
 [ 0.2  0.4]
 [ 0.   0.6]]

Answer 4

One way around it would be to define a different inequality comparison function,

function less_or_nearly_equal(a,b)
    return (a <= b + .000001)

and use that instead of <= in your comparison

for p11,p22 in product(p1,p2):
    if less_or_nearly_equal(p11+p22,max(p1)):    
        P.append((p11,p22))

(I hope my syntax is correct in the function definition.)