Matching Python modulus/remainder in C

Question

Is there a C or C++ function (or compiler flag) for the modulo/modulus/remainder operation (% or %%) for floating point numbers that would conform to the same mathematical standards in corner cases as the equivalent function/operator in Python or R does?

Normally, the operation should be something like

double modulus(x, y)
{
    return x - (y * floor(x / y));
}

But there are many cases in which such simple operation would not yield the same result as Python and R.

These are some examples in Python and R, which do not match with either fmod or remainder or the macro above, even after adding -frounding-math and -fsignaling-nans:

1 % (-1/3) = -1/3
1 % (-1/2) = 0
2 % Inf = 2
(-2) % Inf = Inf
2 % (-Inf) = -Inf
(-2) % (-Inf) = -2

fmod gives instead:

1 % (-1/3) = 0
1 % (-1/2) = 0
2 % (Inf) = 2
-2 % (Inf) = -2
2 % (-Inf) = 2
-2 % (-Inf) = -2

I figured out some additional conditions could be added in order to mimic the results:

double modulus(x, y)
{
    if (isinf(x)) {
        return NAN;
    }
    
    if (isinf(y)) {
        if (isinf(x))
            return NAN;
        else {
            if (x == 0 || isnan(x) || (x > 0) == (y > 0))
                return x;
            else
                return y;
        }
    }

    return x - (y * floor(x / y));
}

But even then it still doesn't work correctly for the cases in which -1 < y < 0, and I don't even know if I'm still missing more corner cases.

How can I get a function that would give the exact same results in C/C++ as the modulus operators from Python and R?

Some of these might be artifacts? Consider: `fmod (1.0, -1.0/3.0) = 5.5511151231257827e-17`, but `fmodf (1.0f, -1.0f/3.0f) = 3.3333331346511841e-01`. — njuffa, Apr 12 '21 at 21:31
@njuffa I guess some of it might be. From the R docs, `mod(x,y) == (x %% y) + y * (x %/% y)`, and `1 %/% (-1/3) = -4`, which is clearly wrong and that's where the mismatch comes. Nevertheless, if both R and Python make the same mistake I guess they are using a common C function. — anymous.asker, Apr 12 '21 at 22:08
`x - (y * floor(x / y))` will very likely be incorrect, if `abs(x)` is much bigger than `abs(y)`, i.e., if `floor(x/y)` can't be encoded as `double`. — chtz, Apr 12 '21 at 22:41
@anymous.asker If the R and Python functionality somehow maps back to the standard C library (which is a reasonable hypothesis), you should be able to lift the relevant code directly from their respective source bases: "Use the source, Luke". — njuffa, Apr 12 '21 at 22:42

score 0 · Answer 1 · answered Apr 12 '21 at 23:07

I didn't manage to find the source code for what Python does, and the NumPy source code (which matches with Python) as far as I understood calls SIMD instructions according to different types, which is hard to read, but I managed to find the R source code for their modulus.

In the end, this doesn't call any built-in C function or assembly instructions, and I'd suspect it could have been taken from Cephes a long time ago before the C99 standard came.

/* Equivalences in order to work with C */
#define R_NaN NAN
#define c_eps DBL_EPSILON
#define R_FINITE(x) (!isinf(x))
typedef long double LDOUBLE;
#define _(msg) (msg)
#define warning(msg) (msg)

/* The actual code */
static double myfmod(double x1, double x2)
{
    if (x2 == 0.0) return R_NaN;
    if(fabs(x2) * c_eps > 1 && R_FINITE(x1) && fabs(x1) <= fabs(x2)) {
    return
        (fabs(x1) == fabs(x2)) ? 0 :
        ((x1 < 0 && x2 > 0) ||
         (x2 < 0 && x1 > 0))
         ? x1+x2  // differing signs
         : x1   ; // "same" signs (incl. 0)
    }
    double q = x1 / x2;
    if(R_FINITE(q) && (fabs(q) * c_eps > 1))
    warning(_("probable complete loss of accuracy in modulus"));
    LDOUBLE tmp = (LDOUBLE)x1 - floor(q) * (LDOUBLE)x2;
    return (double) (tmp - floorl(tmp/x2) * x2);
}

Taken from: https://github.com/wch/r-source/blob/trunk/src/main/arithmetic.c

score -1 · Answer 2 · answered Apr 13 '21 at 14:28

How can I get a function that would give the exact same results in C/C++ as the modulus operators from Python and R?

This is a curious goal as usually the goal is based on an objective standard, not behave just like another language, even if that language does not well perform.

In any case, C has fmod(), a remainder function for OP's consideration.

C math functions are not specified to give a particular answer. sin(x) on one implementation may differ from sin(x) an another, the goal being to match well the math sine(x). How well that is done is a quality of implementation and not a C spec.

A good fmod(x,y) can be expected to be exact. OP's modulus() suffers from imprecision and range issues with x - (y * floor(x / y));

Matching Python modulus/remainder in C

2 Answers2