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

Question

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.

Answer 1

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, http://www.intel.com/design/pentium/manuals/24143004.pdf). 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.

Answer 2

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);
}    

Answer 3

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