CODEWITHSUNDEEP
pythonarraysnumpy
MikeL
MikeL

Reputation: 2469

Equivalent of `math.remainder` in NumPy?

What is the equivalent of math.remainder() function in NumPy? Basically I would like to compute y = x - np.around(x) for a NumPy array x. (It is important that all elements of y have absolute at most 1/2). Looking at NumPy documentation neither np.fmod nor np.remainder do the job.

Obviously, I thought about writing x - np.around(x) but I am afraid that for large values of x subtraction produces floating-point inaccuracies. For example:

import numpy as np

a = np.arange(1000) * 1e-9
x = a / 1e-9
y = x - np.around(x)

Should produce an all-zero vector y, but in practice there will be some errors (that get larger if I increase the size of arrays from 1000 to 10000).

The reason I am asking this question is to figure out if there is a NumPy function for this purpose that calls directly C math library remainder (as math.remainder must do) so as to minimize the floating-points errors.

Upvotes: 2

Views: 498

Answers (2)

aka.nice
aka.nice

Reputation: 9572

IMO, vectorize is the correct solution.
However, if it's considered too slow, I suggest using this:

  • First get the floor modulo
  • then the centered modulo by subtracting around or rint

Like this:

def np_remainder(x,y)
    positive_modulo = x % y
    centered_modulo = positive_modulo - y * numpy.rint(positive_modulo / y)
    return centered_modulo

I think that numpy.fmod is exact, and using it avoid the problems of large rint values.
Since rint will be either 0 or 1 for y>0, the second operation will also be exact (either subtract 0 or y, which should be exact thanks to Sterbenz lemna).
Hence, the formulation above should not introduce any inaccuracy, unlike the raw formulation.

The difference with numpy.vectorize(math.remainder) will be on edge cases of exact tie:

  • The workaround I'm proposing will always return a positive y/2 for y>0
  • math.remainder will choose to round the quotient to nearest even and alternate -y/2 and y/2

Upvotes: 0

Sam Mason
Sam Mason

Reputation: 16213

I don't think this currently exists numpy. That said, the automagic numpy.vectorize call does the right thing for me, e.g:

import math
import numpy as np

ieee_remainder = np.vectorize(math.remainder)

ieee_remainder(np.arange(5), 5)

this broadcasts over the parameters nicely, giving:

array([ 0.,  1.,  2., -2., -1.])

which might be what you want.

performance is about 10 times slower than you'd get with a native implementation. given an array of 10,000 elements my laptop takes approximately the following times:

  • 1200 µs for ieee_remainder
  • 150 µs for a Cython version I hacked together
  • 120 µs for a C program doing a naive loop
  • 80 µs for numpy.fmod

given the complexity of glibc's remainder implementation I'm amazed it's as fast as it is.

Upvotes: 2

Related Questions