How to avoid incorrect rounding with numpy.round?

Question

I'm working with floating point numbers. If I do:

import numpy as np
np.round(100.045, 2)

I get:

Out[15]: 100.04

Obviously, this should be 100.05. I know about the existence of IEEE 754 and that the way that floating point numbers are stored is the cause of this rounding error.

My question is: how can I avoid this error?

Answer 1

You are partly right, often the cause of this "incorrect rounding" is because of the way floating point numbers are stored. Some float literals can be represented exactly as floating point numbers while others cannot.

>>> a = 100.045
>>> a.as_integer_ratio()  # not exact
(7040041011254395, 70368744177664)

>>> a = 0.25
>>> a.as_integer_ratio()  # exact
(1, 4)

It's also important to know that there is no way you can restore the literal you used (100.045) from the resulting floating point number. So the only thing you can do is to use an arbitrary precision data type instead of the literal. For example you could use Fraction or Decimal (just to mention two built-in types).

I mentioned that you cannot restore the literal once it is parsed as float - so you have to input it as string or something else that represents the number exactly and is supported by these data types:

>>> from fractions import Fraction
>>> f = Fraction(100045, 100)
>>> f
Fraction(20009, 20)

>>> f = Fraction("100.045")
>>> f
Fraction(20009, 20)

>>> from decimal import Decimal
>>> Decimal("100.045")
Decimal('100.045')

However these don't work well with NumPy and even if you get it to work at all - it will almost certainly be very slow compared to basic floating point operations.

>>> import numpy as np

>>> a = np.array([Decimal("100.045") for _ in range(1000)])
>>> np.round(a)
AttributeError: 'decimal.Decimal' object has no attribute 'rint'

In the beginning I said that you're are only partly right. There is another twist!

You mentioned that rounding 100.045 will obviously give 100.05. But that's not obvious at all, in your case it is even wrong (in the context of floating point math in programming - it would be true for "normal calculations"). In many programming languages a "half" value (where the number after the decimal you're rounding is 5) isn't always rounded up - for example Python (and NumPy) use a "round half to even" approach because it's less biased. For example 0.5 will be rounded to 0 while 1.5 will be rounded to 2.

So even if 100.045 could be represented exactly as float - it would still round to 100.04 because of that rounding rule!

>>> round(Fraction("100.045"), 1)
Fraction(5002, 5)

>>> 5002 / 5
1000.4

>>> d = Decimal("100.045")
>>> round(d, 2)
Decimal('100.04')

This is even mentioned in the NumPy docs for numpy.around:

Notes

For values exactly halfway between rounded decimal values, NumPy rounds to the nearest even value. Thus 1.5 and 2.5 round to 2.0, -0.5 and 0.5 round to 0.0, etc. Results may also be surprising due to the inexact representation of decimal fractions in the IEEE floating point standard [R1011] and errors introduced when scaling by powers of ten.

(Emphasis mine.)

The only (at least that I know) numeric type in Python that allows setting the rounding rule manually is Decimal - via ROUND_HALF_UP:

>>> from decimal import Decimal, getcontext, ROUND_HALF_UP
>>> dc = getcontext()
>>> dc.rounding = ROUND_HALF_UP
>>> d = Decimal("100.045")
>>> round(d, 2)
Decimal('100.05')

Summary

So to avoid the "error" you have to:

• Prevent Python from parsing it as floating point value and
• use a data type that can represent it exactly
• then you have to manually override the default rounding mode so that you will get rounding up for "halves".
• (abandon NumPy because it doesn't have arbitrary precision data types)

Answer 2

Basically there is no general solution for this problem IMO, unless you have a general rule for all the different cases (see Floating Point Arithmetic: Issues and Limitation). However, in this case you can round the decimal part separately:

In [24]: dec, integ = np.modf(100.045)

In [25]: integ + np.round(dec, 2)
Out[25]: 100.05

The reason for such behavior is not because separating integer from decimal part makes any difference on round()'s logic. It's because when you use fmod it gives you a more realistic version of the decimal part of the number which is actually a rounded representation.

In this case here is what dec is:

In [30]: dec
Out[30]: 0.045000000000001705

And you can check that round gives same result with 0.045:

In [31]: round(0.045, 2)
Out[31]: 0.04

Now if you try with another number like 100.0333, the decimal part is a slightly smaller version which as I mentioned, the result you want depends on your rounding policies.

In [37]: dec, i = np.modf(100.0333)

In [38]: dec
Out[38]: 0.033299999999997

There are also modules like fractions and decimal that provide support for fast correctly-rounded decimal floating point and rational arithmetic, that you can use in situations as such.