How can I simulate GNU C Library drem / remainder function in Google App Engine Python 2.7 runtime?

Question

The GNU C Library has the function drem (alias remainder).

How can I simulate this function just using the modules supported by Google App Engine Python 2.7 runtime?

From the GNU manual for drem:

These functions are like fmod except that they round the internal quotient n to the nearest integer instead of towards zero to an integer. For example, drem (6.5, 2.3) returns -0.4, which is 6.5 minus 6.9.

From the GNU manual for fmod:

These functions compute the remainder from the division of numerator by denominator. Specifically, the return value is numerator - n * denominator, where n is the quotient of numerator divided by denominator, rounded towards zero to an integer. Thus, fmod (6.5, 2.3) returns 1.9, which is 6.5 minus 4.6.

Reading the documentation the following Python code should work:

def drem(x, y):
    n = round(x / y)
    return x - n * y

However with Python, drem(1.0, 2.0) == -1.0 and with C drem(1.0, 2.0) == 1.0. Note Python returns negative one and C returns positive one. This is almost certainly an internal difference in rounding floats. As far as I can tell both functions perform the same otherwise where parameters 2 * x != y.

How can I make my Python drem function work the same as its C equivalent?

Answer 1

The key to solving this problem is to realise that the drem/remainder function specification requires the internal rounding calculation to round to half even.

Therefore we cannot use the built-in round function in Python 2.x as it rounds away from 0. However the round function in Python 3.x has changed to round to half even. So the following Python 3.x code will be equivalent to the GNU C Library drem function but will not work in Python 2.x:

def drem(x, y):
    n = round(x / y)
    return x - n * y

To achieve the same with Python 2.x we can use the decimal module and its remainder_near function:

import decimal

def drem(x, y):
    xd = decimal.Decimal(x)
    yd = decimal.Decimal(y)

    return float(xd.remainder_near(yd))

Answer 2

EDIT: I just read your first comment and see that you cannot use the ctypes module. Anyways, I learned a lot today by trying to find an answer to your problem.

Considering that numpy.round() rounds values exactly halfway between rounded decimal values to the next even integer, using numpy is not a good solution.

Also, drem internally calls this MONSTER function, which should be hard to implement in Python.

Inspired by this article, I would recommend you to call the drem function from the math library directly. Something along these lines should do the trick:

from ctypes import CDLL
# Use the C math library directly from Python
# This works for Linux, but the version might differ for some systems
libm = CDLL('libm.so.6') 

# For Windows, try this instead: 
# from ctypes import cdll
# libm = cdll.libc

# Make sure the return value is handled as double instead of the default int
libm.drem.restype = c_double
# Make sure the arguments are double by putting them inside c_double()
# Call your function and have fun!
print libm.drem(c_double(1.0), c_double(2.0))