Modulo function in c++, that behaves like mod in matlab

I do a simulation with a lot of particles (up to 100000) in periodic domain(box), and in order particles to stay inside the box, I use modulo function with float or double numbers.

In Matlab everything works great with mod function. However in C++ I found out, that function fmod is not completely equal to Matlab's mod function:

mod(-0.5,10)=9.5 - I want this result in C++

fmod(-0.5,10)=-0.5 - I don't want this.

I, of course, can solve my problem with if statements. However, i think, it will affect efficiency (if statement in critical loop). Is there a way to implement this function without if statement? May be some other function?

Thanks.

5 Answers

Just use a conditional. It will not meaningfully affect efficiency.

inline double realmod (x, y) { result = fmod(x, y); return result >= 0 ? result : result + y; } 

fmod() calls assembly instruction FPREM which takes 16-64 cycles (according to the Pentium manual, ). The jump instructions for the conditional and the floating point addition only amount to 5 or so.

When your code has floating point division, you don't need to sweat the small stuff.

Either use floor and regular division:

float modulo(float a, float q) { float b = a / q; return (b - floor(b)) * q; } 

or you can add the divisor to the result of fmod without branching:

float modulo(float a, float q) { float m = fmod(a, q); return m + q * (m < 0.f); } 

Based on Matlab mod(a, m) documentation and @QuestionC's answer -

A general solution that behaves exactly like Matlab - also for negative and zero divisor.

Tested against multiple values :

static inline double MatlabMod(double q, double m) { if(m == 0) return q; double result = fmod(q, m); return ((result >= 0 && m > 0) || (q <= 0 && m < 0)) ? result : (result + m); } 

Tested with matlab for :

(q, m) -> result

(54, 321) -> 54

(-50, 512) -> 462

(54, -152) -> -98

(-53, -500) -> -53

(-500, 300) -> 100

(-5000, 400) -> 200

(-1000, -360) -> -280

(500, 360) -> 140

(1000, 360) -> 280

(-1000, 360) -> 80

(-5051, 0) -> -5051

(512, 0) -> 512

(0, 52) -> 0

(0, -58) -> 0

Just add the divisor to the number you want to keep in the interval before you apply the modulo operator:

return fmod(a+q,q); 

this requires no branching at all.

If you have to worry about a exeeding -q between two updates, you can make it more robust by e.g.:

return fmod(a+q*10,q); 

which would work for a >= -10*q

The most straightforward with working for both floats and ints without any branching.

// b = MOD(a, m) int a = a - m * floor(a / m) + m; int b = a - m * floor(a / m);