Why does python round(np.float16(np.pi),5) return infinity? Bug, limitation, or expected?

Question

I have a Python 3 function rounding floats to 6 digits (the logic handles various precision levels). When passed with many (possibly all) numpy.float16 values it generates multiply-overflow warnings and returns infinity.

The short snippet in the question title or shown below illustrates the behavior.

The workaround is easy, just convert to larger floats first, but I'm curious whether the behavior is expected.

import numpy as np
x = np.float16(3.14)
x = round(x, 5)
if np.isinf(x):
    print("you've made an infinity through rounding....", 1, x)
else:
    print('just x: ', x)

I expect round to affect precision, but never cause an overflow or change a value into infinity.

Answer 1

This is a limitation of NumPy's round algorithm. I hesitate to call it a bug: that's for the NumPy core developers to decide, but it may be worth reporting nevertheless.

Here's the problem: to round to 5 decimal places, NumPy does the equivalent of scaling by 100000.0, rounding to the nearest integer, then diving by 100000.0 again. That initial scaling can potentially overflow, even when the final result of the round operation would be expected to be within range.

Here's the portion of the NumPy source where this is implemented. You need to backtrack a bit in the source to figure out that in this case op1 is multiplication and op2 refers to division.

With float64 or float32 this would be unlikely to cause issues, because for normal uses you're unlikely to be within 100000.0 of the upper bound of the representable range of the floating-point type. But if you do get too close to that upper bound, you'll see the same problem. Here's an example with np.float64:

>>> np.finfo(np.float64).max
1.7976931348623157e+308
>>> x = np.float64(1e304)  # pick something within 1e5 of that max
>>> x * 1e5
__main__:1: RuntimeWarning: overflow encountered in double_scalars
inf
>>> np.round(x, 5)  # same multiplication happening internally
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/numpy/core/fromnumeric.py:56: RuntimeWarning: overflow encountered in multiply
  return getattr(obj, method)(*args, **kwds)
inf

And here's the same thing with float32:

>>> np.finfo(np.float32).max
3.4028235e+38
>>> x = np.float32(1e35)
>>> x * 1e5  # okay; NumPy converts to `float64`
1.0000000409184788e+40
>>> np.round(x, 5)
inf

With np.float16, it's exactly the same issue, but since the dynamic range of the float16 type is so small, you're much more likely to observe the issue in practice.

In general, though, note that even if this were fixed, it is possible for two-argument round to overflow: it's possible that the original value is within the range of the relevant floating-point type, while the rounded value is not. Here's an example with Python's own round function:

>>> x = 1.76e308
>>> x
1.76e+308
>>> round(x, -307)  # should be 1.8e308, but that's out of range
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
OverflowError: rounded value too large to represent

But this can only happen for a negative ndigits argument. If the second argument were nonnegative, then there's no possibility of overflow - every sufficiently large representable value in any of the standard floating-point types is already integral, so a round with a nonnegative ndigits shouldn't change its value.