http://ridiculousfish.com/blog/posts/labor-of-division-episo...

http://ridiculousfish.com/blog/posts/labor-of-division-episo...

http://ridiculousfish.com/blog/posts/labor-of-division-episo...

The same author wrote a library for doing this (and more; it also uses bit shifts) at runtime: http://libdivide.com

Of course, that only makes sense if you know you will do lots of divisions by the same number

http://www.hackersdelight.org/

https://github.com/dlang/dmd/blob/master/src/backend/divcoef...

[1] https://doi.org/10.1145/773473.178249

LLVM does not currently generate that code, but you can make a case it should.

    λ> import Data.SBV
    λ> prove $ do x <- Data.SBV.sWord32 "x"; return $ (x * 0xaaaaaaab .< 0x55555556) .== (x `sMod` 3 .== 0)
    Q.E.D.

    λ> sat $ do [a,b] <- sWord16s ["a", "b"]; forAll_ (\x -> (x * a .< b) <=> (x `sMod` 7 .== 0))
    Satisfiable. Model:
      a = 28087 :: Word16
      b =  9363 :: Word16

    λ> sat $ do [a,b] <- sWord16s ["a", "b"]; forAll_ (\x -> (x * a .< b) <=> (x `sMod` 9 .== 0))
    Satisfiable. Model:
      a = 36409 :: Word16
      b =  7282 :: Word16

    λ> sat $ do [a,b] <- sWord16s ["a", "b"]; forAll_ (\x -> (x * a .< b) <=> (x `sMod` 11 .== 0))
    Satisfiable. Model:
      a = 35747 :: Word16
      b =  5958 :: Word16

Every uint32 can be expressed as k=3*N mod 2^32. We can get that N by multiplying by the modular inverse. k is divisible by 3 if and only if the multiplication 3N doesn't "wrap around".

So: Multiply by modular inverse to get N, check that it's small enough that 3N < 2^32 ore equivalently that N < 2^32 / 3 = 0x55555556.

    gcd15(int x) {
         return ((x * 0xaaaaaaab) < (0xfffffff/3 + 1) ? 1 : 3)
              * ((x * 0xaaaaaaab) < (0xfffffff/5 + 1) ? 1 : 5)
    }

    gcd15(int x) { 
       LOOKUP_TABLE = { /* gcd(n,15) for n = 0,...,15 */ }
       int q = (x * ((1 << 28) / 15)) >> 28;
       return LOOKUP_TABLE[x - 15 * q];
    }

Disclaimer: That rule is restricted to fast processors, not low-power embedded stuff.

Inlining is not that important for speed itself. It rather is an optimization booster, which makes many other optimization much more effective. If your other optimizations are crap, then inlining will not help much, because the call overhead is not that big.

There's a nice overview for people with a passing interest here: http://www.lighterra.com/papers/basicinstructionscheduling/

I might start prototyping some things in DAG notation just to avoid ever forgetting this conceptual abstraction.

(I'm one of the authors)

Why can't one just multiply with the inverse of 19 (which can be calc'ed during compile time)?

1. Bit shifting

2. Integer multiplication

3. Integer division

4. Floating point multiplication

The trick in the article works because the cost of 1 + 2 is still smaller than 3.

Multiplying with the inverse of 19, a floating point number, wouldn't work because 4 is more costly than 3.

[0] http://nicolas.limare.net/pro/notes/2014/12/12_arit_speed/

> IMPORTANT: Useful feedback revealed that some of these measures are seriously flawed. A major update is on the way.

Looking over the results, some of the numbers are off.

On Intel CPUs, FP multiplication is faster than integer division. Might not be true on ARM CPUs which generally have slower FPUs.

On Skylake, for example, 32-bit unsigned integer division has a 26 cycle latency with a throughput of 1 instruction / 6 cycles, while 32/64-bit floating point multiplication has a 4 cycle latency with a throughput of 2 instructions / cycle.

Source: http://agner.org/optimize/instruction_tables.pdf

For divisions by a constant value that don't easily decompose into shifts you can fall back to multiplication by a magic constant which is the integer reciprocal. (This is also something compilers do and is what's being explained in the article.)

What's being explained in the article is multiplying by a fraction the value of which is close to the rational reciprocal of the divisor, and where the denominator of the fraction is an integer power of two (so dividing by the denominator can be done with a shift).

The fraction in this case is (2938661835 + 2^32) / 2^37.

The approximate latencies for Skylake are:

    div --> 26 cycles
    cvtsi2sd + mulsd + cvttsd2siq --> 6 + 4 + 6 = 16 cycles

    Real integer division (-Os) --> 1.392s
    FPU Multiply (-Os)          --> 0.243s
    FPU Multiply (-O2)          --> 0.197s
    Integer Multiply (-O2)      --> 0.164s

    #include <stdio.h>

    int main() {
        volatile unsigned x;
        for (unsigned n = 0; n < 100000000; ++n) {
    #if 1 /* Change to 0 to use FPU. */
            /*
            Compile with -Os to get GCC to emit div instruction.
            -O2 to emit integer multiply.
            Clang emits integer multiply, even with -Os.
            */
            x = n / 19;
    #else
            /* Use the FPU. */
            x = (double)n * (1.0 / 19.0);
    #endif
        }
    }

For n < 19, "(double)n * (1.0 / 19.0)" evaluates to a double between 0.0 and 1.0, then it is truncated to 0 when it is implicitly converted to unsigned int.

Since there are only 2^32 values for 32-bit integers, it is possible to test all values in under a minute:

    #include <stdio.h>
    #include <stdint.h>

    int main() {
        uint32_t n = 0;
        do {
            uint32_t a = n / 19;
            uint32_t b = (double)n * (1.0 / 19.0);
            if (a != b) {
                printf("Not equal for n = %u\n", n);
            }
            ++n;
        } while (n != 0);
    }

If you round the FP division up (to the next largest representable value, that is), it should be correct, for 32 bit integer types at least.

> On Skylake [...] 32/64-bit floating point multiplication has a 4 cycle latency with a throughput of 2 instructions / cycle.

Of course, there are some operations that are very expensive (trigonometric functions, for example), but they're not necessary here, and they're also very expensive on the GPU.

[1] https://news.ycombinator.com/item?id=13693749

This is wrong, the divisor and dividend are the wrong way around.

..and in case anyone else is wondering how those numbers relate to the constant in the code, it's

  2^37 / (2938661835 + 2^32)

19 * 678,152,731 = 1 (mod 2^32)

19 * 9,708,812,670,373,448,219 = 1 (mod 2^64)

As long as we already know the size of an int Euclid's algorithm should be able to find these inverses if they exist.

What algorithm do you have in mind?

[1] https://en.wikipedia.org/wiki/Pentium_FDIV_bug