BUG rounding error in divmod

The ndarray implementation of divmod seems to do some incorrect rounding:

# regular python
divmod(78*6e-8, 6e-8)
# (77.0, 5.999999999999965e-08)
# makes sense; cannot represent number precisely as float
import numpy as np
np.__version__
# '1.10.0.dev0+00f4fae'
divmod(np.array(78*6e-8), 6e-8)
# (78.0, 5.9999999999999651e-08)
# Oops!!
divmod(np.arange(77, 80)*6e-8, 6e-8)
# (array([ 77.,  78.,  79.]),
#   array([  2.24993127e-22,   6.00000000e-08,   6.00000000e-08]))

@MatthewQuenneville

p.s. Bug initially found by @MatthewQuenneville

A bit further checking shows that divmod just calls floor_divide and remainder separately
(https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/number.c#L781), and that floor_divide gives a different answer from python's version:

78*6e-8 // 6e-8
# 77.0
np.floor_divide(78*6e-8, 6e-8)
# 78.0

Of course, what is "wrong" is a bit in the eye of the beholder:

78*6e-8 / 6e-8 == 78
# True

In any case, at least the divmod should give a remainder that is consistent with the fraction. I fear one either has to write a new inner loop or ensure the remainder is calculated using the fraction (or vice versa).

A possible cause of the problem may be a check on negative remainders done in https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/loops.c.src#L1004: the remainder is increased by one unit, but the fraction would not be correspondingly decreased.

Given that:

>>> (78*6e-8) / 6e-8
78.0

It seems that Python's floor division is a bit problematic:

>>> (78*6e-8) // 6e-8
77.0

For reference, I reported a Python issue at http://bugs.python.org/issue25129

Not really, that could simply be a floating point rounding error, and the remainder is consistent with that.

All that matters is that one can add back, and this is fine for python:

q, r = divmod(78*6e-8, 6e-8)
q, r
# (77.0, 5.999999999999965e-08)
(q * 6e-8 + r) / 6e-8
# 78.0

In contrast, with a numpy array:

q, r = divmod(np.array(78*6e-8), 6e-8)
(q * 6e-8 + r) / 6e-8
# OOPS:  78.999999999999986

Looks to me like it depends on the details of the computation

In [6]: '%25.18e' % ((78 * 6e-8)/6e-8)
Out[6]: ' 7.800000000000000000e+01'

In [7]:  a = (78 * 6e-8)

In [8]:  a / 6e8
Out[8]: 7.799999999999999e-15

nvm. Yes, it looks like numpy gets the remainder wrong.

nvm. Yes, it looks like numpy gets the remainder wrong.

I think the problem is that the two are not calculated together, so that the adjustment in the remainder is not propagated to the quotient.

Indeed, CPython internally has a unique divmod function that gets called by all three operations: divmod(), floor division and remainder.

Consider a hypothetical floor_sub(78, 1e-15). Would you want it to give 78 because 78 - 1e-15 == 78 in double precision or would you want it to give the more mathematically correct answer 77?

@argriffing

@argriffing - I know one cannot always avoid precision/rounding errors; however, the remainder and floor division should be consistent with each other, in that if q, r = divmod(a, b), one should have q*b+r = a, at least approximately; this is violated in the numpy example I showed.

@mhvk

@mhvk Yes, I was responding to @pitrou's suggestion that Python's floor division is problematic. Sorry that wasn't clear in my comment! I agree that numpy has a divmod bug.

Here's a case that's perhaps less contrieved, only using "non-small" numbers:

In [3]: divmod(np.float64(1.0), 0.2)
Out[3]: (5.0, 0.19999999999999996) # oops

I think we can get around the same path problem as long as division and floor are repeatable. Which may not always be the case on x86, but...

Happy to see this resolved; solution looks good!

@mhvk

@mhvk, I think it is actually better than Python at the moment. However, at some point it will hit the precision limits of the floats...

There seems to be a regression now: #7224

More precise choice of axes limits.

BUG: Make divmod behave better under roundoff error.

TST: skip broken test

Navigation Menu

Repository navigation

BUG rounding error in divmod #6127

20 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

BUG rounding error in divmod #6127

Description

Activity

mhvk commented on Jul 28, 2015

mhvk commented on Jul 29, 2015

pitrou commented on Sep 15, 2015

pitrou commented on Sep 15, 2015

mhvk commented on Sep 15, 2015

charris commented on Sep 15, 2015

charris commented on Sep 15, 2015

mhvk commented on Sep 15, 2015

pitrou commented on Sep 15, 2015

argriffing commented on Sep 15, 2015

mhvk commented on Sep 15, 2015

argriffing commented on Sep 15, 2015

anntzer commented on Dec 30, 2015

20 remaining items

charris commented on Jan 15, 2016

mhvk commented on Jan 17, 2016

charris commented on Jan 17, 2016

pitrou commented on Feb 10, 2016

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions