Fast Euclidean division in C

Question

I am interested in getting the remainder of the Euclidean division, that is, for a pair of integers (i, n), find r such as:

i = k * n + r, 0 <= r < |k|

the simple solution is:

int euc(int i, int n)
{
    int r;

    r = i % n;
    if ( r < 0) {
        r += n;
    }
    return r;
}

But since I need to execute this tens of million of times (it is used inside an iterator for multidimensional arrays), I would like to avoid the branching if possible. Requirements:

• Branching but faster is also desirable.
• A solution which works only for positive n is acceptable (but it has to work for negative i).
• n is not known in advance, and can be any value > 0 and < MAX_INT

Edit

It is actually quite easy to get the result wrong, so here is an example of the expected results:

• euc(0, 3) = 0
• euc(1, 3) = 1
• euc(2, 3) = 2
• euc(3, 3) = 0
• euc(-1, 3) = 2
• euc(-2, 3) = 1
• euc(-3, 3) = 0

Some people also worry that it does not make sense to optimize this. I need this for an multi-dimensional iterator where out of bounds items are replaced by items in a 'virtual array' which repeats the original array. So if my array x is [1, 2, 3, 4], the virtual array is [...., 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4], and for example, x[-2] is x[1], etc...

For a nd array of dimension d, I need d Euclidean division for every point. If I need to do a correlation between a n^d array with a m^d kernel, I need n^d * m^d * d euclidean divisions. For a 3d image of 100x100x100 points and a kernel of 5*5*5 points, that's already ~ 400 million Euclidean divisions.

Answer 1

Edit: No multiplication or branches woot.

int euc(int i, int n)
{
    int r;

    r = i % n;
    r += n & (-(r < 0));

    return r;
}

Here is the generated code. According to the MSVC++ instrumenting profiler (my testing) and the OP's testing, they perform nearly the same.

; Original post
00401000  cdq              
00401001  idiv        eax,ecx 
00401003  mov         eax,edx 
00401005  test        eax,eax 
00401007  jge         euc+0Bh (40100Bh) 
00401009  add         eax,ecx 
0040100B  ret              

; Mine
00401020  cdq              
00401021  idiv        eax,ecx 
00401023  xor         eax,eax 
00401025  test        edx,edx 
00401027  setl        al   
0040102A  neg         eax  
0040102C  and         eax,ecx 
0040102E  add         eax,edx 
00401030  ret

Answer 2

int euc(int i, int n)
{
    return (i % n) + (((i % n) < 0) * n);
}

Answer 3

"Premature optimization is the root of all evil"

Answer 4

If you can also guarantee that i is never less than -n, you can simply put the optional addition before the modulo. That way, you don't need the branch, and the modulo cuts out what you added if you didn't need to.

int euc(int i, int n)
{
    return (i + n) % n;
}

If i is less than -n, you can still use this method. In a case like this, you probably know exactly what range your values will be in. Therefore, instead of adding n to i, you can add x*n to i where x is any integer that gives you a sufficient range. For added speed (on processors that don't have single-cycle multiply), you could left-shift instead of multiplying.

Answer 5

Integer multiplication is much faster than division. For a large number of calls with a known N, you can replace division by N by multiplication by a pseudo-inverse of N.

I will illustrate this on an example. Take N=29. Then compute once a pseudo inverse 2^16/N: K=2259 (truncated from 2259.86...). I assume I is positive and I*K fits on 32 bits.

Quo = (I*K)>>16;   // replaces the division, Quo <= I/N
Mod = I - Quo*N;   // Mod >= I%N
while (Mod >= N) Mod -= N;  // compensate for the approximation

In my example, let's take I=753, we get Quo=25 and Mod=28. (no compensation needed)

EDIT.

In your 3D convolution example, most of the calls to i%n will be with i in 0..n-1, so in most cases, a first line like

if (i>=0 && i<n) return i;

will bypass the costly and here useless idiv.

Also, if you have enough RAM, just align all dimensions to powers of 2 and use bit manipulations (shift,and) instead of divisions.

EDIT 2.

I actually tried it on 10^9 calls. i%n: 2.93s, my code: 1.38s. Just keep in mind it implies a bound on I (I*K must fit on 32 bits).

One other thought: if your values are x+dx, with x in 0..n-1 and dx small, then the following will cover all the cases:

if (i<0) return i+n; else if (i>=n) return i-n;
return i;

Answer 6

I timed everyone's proposals in gcc -O3 using TSC (except the one for constant N), and they all took the same amount of time (within 1%).

My thought was that either ((i%n)+n)%n (no branching), or (i+(n<<16))%n (obviously fails for large n or extremely negative i) would be faster, but they all took the same time.

Answer 7

For tens of millions this is not worth thinking about. Maybe for tens of billions, but then do proper benchmarking.

Answer 8

I really like the expression:

r = ((i%n)+n)%n;

The disassembly is very short:

r = ((i%n)+n)%n;

004135AC  mov         eax,dword ptr [i] 
004135AF  cdq              
004135B0  idiv        eax,dword ptr [n] 
004135B3  add         edx,dword ptr [n] 
004135B6  mov         eax,edx 
004135B8  cdq              
004135B9  idiv        eax,dword ptr [n] 
004135BC  mov         dword ptr [r],edx

It has no jumps (2 idivs, which might be costly), and it can be completely inlined, avoiding the overhead of a function call.

What do you think?

Answer 9

If you have low enough range create a look-up table - two dim array. Also you can make the function Inline, and make sure it is by looking at the produced code.

Answer 10

I think 280Z28 and Christopher have the assembler golf covered better than I would, and that deals with random access.

What you're actually doing, though, seems to be processing entire arrays. Obviously for reasons of memory caching you already want to be doing this in order if possible, since avoiding a cache miss is a many, many times better optimisation than avoiding a small branch.

In that case, with a suitable bounds check first you can do the inner loop in what I will call "dashes". Check that the next k increments don't result in an overflow in the smallest dimension on either array, and then "dash" k steps using a new even-more-inner loop which just increments the "physical" index by 1 each time instead of doing another idiv. You or the compiler can unroll this loop, use Duff's device, etc.

If the kernel is small, and especially if it is fixed size, then that (or a multiple of it with suitable unrolling to occasionally subtract instead of adding), is probably the value to use for the length of the "dash". Compile-time constant dash length is probably best, since then you (or the compiler) can fully unroll the dash loop and elide the continuation condition. As long as this doesn't make the code too big to be fast, it essentially replaces the entire positive-modulo operation with an integer increment.

If the kernel is not fixed size, but is often very small in its last dimension, consider having different versions of the comparison function for the most common sizes, with the dash loop fully unrolled in each.

Another possibility is to calculate the next point at which an overflow will occur (in either array), and then dash to that value. You still have a continuation condition in the dash loop, but it goes as long as possible using only increments.

Alternatively, if operation you're doing is numeric equality or some other simple operation (I don't know what a "correlation" is) you could look at SIMD instructions or whatever, in which case the dash length should be a multiple of the widest single-instruction compare (or appropriate SIMD op) on your architecture. This isn't something I have experience with, though.

Answer 11

Without a branch, but a bit bit-fiddling :

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (((unsigned int)r) >> 31) * n;
    return r;
}

Without multiply :

int euc2(int i, int n)
{
    int r;
    r = i % n;
    r += (r >> 31) & n;
    return r;
}

This gives :

; _i$ = eax
; _n$ = ecx

cdq
idiv   ecx
mov eax, edx
sar eax, 31
and eax, ecx
add eax, edx

Answer 12

You say in your answer to Eric Bainville that most of the time 0 <= i < n and that you have

if (i>=0 && i<n) return i;

as the first line of your euc() anyway.

As you are making the comparisons anyway, you might as well use them:

int euc(int i, int n)
{
    if (n <= i)            return i % n;
    else if (i < 0)        return ((i + 1) % n) + n - 1;
    else /* 0 <= i < n */  return i;  // fastest possible response for common case
}

Answer 13

If you can guarantee that the dimensions of your array are always powers of two, then you can do this:

r = (i & (n - 1));

If you can further guarantee that your dimensions will be from a given subset, you can do:

template<int n>
int euc(int i) {
    return (i & (n - 1));
}

int euc(int i, int n) {
    switch (n) {
        case 2: return euc<2>(i);
        case 4: return euc<4>(i);
    }
}

Answer 14

Here's Christopher's version with a fallback to Jason's if right shifting is not arithmetic.

#include <limits.h>
static inline int euc(int i, int n)
{
    // check for arithmetic shift
    #if (-1 >> 1) == -1
        #define OFFSET ((i % n >> (sizeof(int) * CHAR_BIT - 1)) & n)
    #else
        #define OFFSET ((i % n < 0) * n)
    #endif

    return i % n + OFFSET;
}

The fallback version should be slower as it uses imul instead of and.