Faster divisibility test than % operator?

Question

I noticed a curious thing on my computer.^* The handwritten divisibility test is significantly faster than the % operator. Consider the minimal example:

^{*_{AMD Ryzen Threadripper 2990WX, GCC 9.2.0}}

static int divisible_ui_p(unsigned int m, unsigned int a)
{
    if (m <= a) {
        if (m == a) {
            return 1;
        }

        return 0;
    }

    m += a;

    m >>= __builtin_ctz(m);

    return divisible_ui_p(m, a);
}

The example is limited by odd a and m > 0. However, it can be easily generalized to all a and m. The code just converts the division to a series of additions.

Now consider the test program compiled with -std=c99 -march=native -O3:

    for (unsigned int a = 1; a < 100000; a += 2) {
        for (unsigned int m = 1; m < 100000; m += 1) {
#if 1
            volatile int r = divisible_ui_p(m, a);
#else
            volatile int r = (m % a == 0);
#endif
        }
    }

... and the results on my computer:

| implementation     | time [secs] |
|--------------------|-------------|
| divisible_ui_p     |    8.52user |
| builtin % operator |   17.61user |

Therefore more than 2 times faster.

The question: Can you tell me how the code behaves on your machine? Is it missed optimization opportunity in GCC? Can you do this test even faster?

UPDATE: As requested, here is a minimal reproducible example:

#include <assert.h>

static int divisible_ui_p(unsigned int m, unsigned int a)
{
    if (m <= a) {
        if (m == a) {
            return 1;
        }

        return 0;
    }

    m += a;

    m >>= __builtin_ctz(m);

    return divisible_ui_p(m, a);
}

int main()
{
    for (unsigned int a = 1; a < 100000; a += 2) {
        for (unsigned int m = 1; m < 100000; m += 1) {
            assert(divisible_ui_p(m, a) == (m % a == 0));
#if 1
            volatile int r = divisible_ui_p(m, a);
#else
            volatile int r = (m % a == 0);
#endif
        }
    }

    return 0;
}

compiled with gcc -std=c99 -march=native -O3 -DNDEBUG on AMD Ryzen Threadripper 2990WX with

gcc --version
gcc (Gentoo 9.2.0-r2 p3) 9.2.0

UPDATE2: As requested, the version that can handle any a and m (if you also want to avoid integer overflow, the test has to be implemented with integer type twice as long as the input integers):

int divisible_ui_p(unsigned int m, unsigned int a)
{
#if 1
    /* handles even a */
    int alpha = __builtin_ctz(a);

    if (alpha) {
        if (__builtin_ctz(m) < alpha) {
            return 0;
        }

        a >>= alpha;
    }
#endif

    while (m > a) {
        m += a;
        m >>= __builtin_ctz(m);
    }

    if (m == a) {
        return 1;
    }

#if 1
    /* ensures that 0 is divisible by anything */
    if (m == 0) {
        return 1;
    }
#endif

    return 0;
}

Comments are not for extended discussion; this conversation has been [moved to chat](https://chat.stackoverflow.com/rooms/209890/discussion-on-question-by-dabler-faster-divisibility-test-than-operator). — Samuel Liew, Mar 19 '20 at 01:12
I would like to also see a test where you actually assert that those two `r`s that you calculate are indeed equal to each other. — Mike Nakis, Mar 19 '20 at 14:28
Most real-life uses of `a % b` have `b` much smaller than `a`. Through most iterations in your test case, they are of similar size, or `b` is bigger, and your version can be faster on many CPUs in those situations. — Matt Timmermans, Mar 19 '20 at 14:37

score 12 · Answer 1 · answered Mar 19 '20 at 14:25

12

What you’re doing is called strength reduction: replacing an expensive operation with a series of cheap ones.

The mod instruction on many CPUs is slow, because it historically was not tested in several common benchmarks and the designers therefore optimized other instructions instead. This algorithm will perform worse if it has to do many iterations, and % will perform better on a CPU where it needs only two clock cycles.

Finally, be aware that there are many shortcuts to take the remainder of division by specific constants. (Although compilers will generally take care of this for you.)

answered Mar 19 '20 at 14:25

Davislor

14,674
2
34
49

*historically was not tested in several common benchmarks* - Also because division is inherently iterative and hard to make fast! x86 at least does remainder as part of `div` / `idiv` which have gotten some love in Intel Penryn, Broadwell and IceLake (higher radix hardware dividers) – Peter Cordes Apr 07 '20 at 11:01
1

My understanding of "strength reduction" is that you replace a heavy operation in a loop with a single lighter operation, e.g. instead of `x = i * const` every iteration you do `x += const` every iteration. I don't think replacing a single multiply with a shift/add loop would be called strength-reduction. https://en.wikipedia.org/wiki/Strength_reduction#Other_strength_reduction_operations says the term can maybe be used this way, but with a note "This material is disputed. It is better described as peephole optimization and instruction assignment." – Peter Cordes Apr 07 '20 at 11:08

score 9 · Accepted Answer · answered Mar 21 '20 at 10:09

I will answer my question myself. It seems that I became a victim of branch prediction. The mutual size of the operands does not seem to matter, only their order.

Consider the following implementation

int divisible_ui_p(unsigned int m, unsigned int a)
{
    while (m > a) {
        m += a;
        m >>= __builtin_ctz(m);
    }

    if (m == a) {
        return 1;
    }

    return 0;
}

and the arrays

unsigned int A[100000/2];
unsigned int M[100000-1];

for (unsigned int a = 1; a < 100000; a += 2) {
    A[a/2] = a;
}
for (unsigned int m = 1; m < 100000; m += 1) {
    M[m-1] = m;
}

which are / are not shuffled using the shuffle function.

Without shuffling, the results are still

| implementation     | time [secs] |
|--------------------|-------------|
| divisible_ui_p     |    8.56user |
| builtin % operator |   17.59user |

However, once I shuffle these arrays, the results are different

| implementation     | time [secs] |
|--------------------|-------------|
| divisible_ui_p     |   31.34user |
| builtin % operator |   17.53user |

Faster divisibility test than % operator?

2 Answers2