Equivalent matlab function mod in numpy or python

Question

I was converting some codes from Matlab to Python.

In matlab there is the function mod which gives the modulo operation.

For example the following example shows different results between the matlab mod and the equivalent numpy remainder operation:

Matlab:

>> mod(6, 0.05)

ans =

     0

Numpy:

np.remainder(6, 0.05) 
0.04999999999999967

np.mod(6, 0.05)
0.04999999999999967

Python modulus operator gives the same results as numpy.

6%0.05
0.04999999999999967

Is there anything in python which gives the same mod operation as the one in Matlab? preferably can be operated on numpy 2d/3d arrays.

numpy documentation says that numpy.mod is equivalent to matlab mod.

Answer 1

This is the core of the problem, in python:

>>> 6/0.05 == 120
True

>>> 6//0.05 == 120  # this is 119 instead
False

The floating-point result of 6/0.05 is close enough to 120 (i.e. within the resolution of double precision) that it gets rounded to 120.0. However, it's ever so slightly smaller than 120, so explicit floor division will truncate that number to 119 before it could be normalized to 120.0.

Some proof:

>>> from decimal import Decimal
... print(6/Decimal(0.05))  # exactly approximate
... print(6/Decimal('0.05'))  # exact
119.9999999999999933386618522
1.2E+2

The first number is what you'd first get with 6/0.05, but the number 119.9999999999999933386618522 gets rounded to the nearest number representable with double precision, and this is 120. One can easily prove that these two numbers are indeed the same within double precision:

>>> print(6/Decimal('0.05') - 6/Decimal(0.05))
6.6613381478E-15

>>> 120 - 6.6613381478E-15 == 120
True

Now here's help mod from MATLAB:

    MOD(x,y) returns x - floor(x./y).*y if y ~= 0, carefully computed to
    avoid rounding error. If y is not an integer and the quotient x./y is
    within roundoff error of an integer, then n is that integer.

This suggests that when x/y is close to an integer then it's rounded first, rather than being truncated like in python. So MATLAB goes out of its way to do some magic with the floating-point results.

The simplest solution would be to round your numbers yourself (unless you can use something like decimal.Decimal, but this means you should forgo native doubles entirely, including literals) and reproduce MATLAB's mod that way, assuming that makes sense for your use cases.

Answer 2

Here is a workaround. It basically rounds the denominator to its str representation and from there does integer arithmetic:

import numpy as np
import decimal

def round_divmod(b,a):
     n,d = np.frompyfunc(lambda x:decimal.Decimal(x).as_integer_ratio(),1,2)(a.astype('U'))
     n,d = n.astype(int),d.astype(int)
     q,r = np.divmod(b*d,n)
     return q,r/d

a = np.round(np.linspace(0.05,1,20),2).reshape(4,5)
a
# array([[0.05, 0.1 , 0.15, 0.2 , 0.25],
#        [0.3 , 0.35, 0.4 , 0.45, 0.5 ],
#        [0.55, 0.6 , 0.65, 0.7 , 0.75],
#        [0.8 , 0.85, 0.9 , 0.95, 1.  ]])
round_divmod(6,a)
# (array([[120,  60,  40,  30,  24],
#        [ 20,  17,  15,  13,  12],
#        [ 10,  10,   9,   8,   8],
#        [  7,   7,   6,   6,   6]]), array([[0.  , 0.  , 0.  , 0.  , 0.  ],
#        [0.  , 0.05, 0.  , 0.15, 0.  ],
#        [0.5 , 0.  , 0.15, 0.4 , 0.  ],
#        [0.4 , 0.05, 0.6 , 0.3 , 0.  ]]))