Item — HN Reader

Here's a trick I came across recently which I found quite neat: you can test for divisibility by 3 using (x * 0xaaaaaaab) < 0x55555556. Same concept works for 5 and 15, but sadly not other factors.

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

3 comments

That's pretty cool, although it ought to be mentioned that x must be a uint32. Here's a proof:

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

Are you sure it doesn't work for other factors? It seems it works for 7, 9, and 11.

    λ> 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

So a is equal to 7^{-1} mod 2^16, and b to ceil((2^16-1)/7), etc.

Duh, you're right on both counts. It should work for any odd divisor, I was just computing the inverse in a dumb way when I tried out other divisors. And I should have stated u32 wrapping arithmetic.

And for even divisors you can separate the divisor into a power of two, and an odd part. The odd part is checked as before, and the power of two is checked with an and mask.

An even trickier way to do this is to rotate right by the power of two factor between doing the multiply and the compare.

That wouldn't work, because a number that is bigger than the threshold could become smaller than it after the rotation.

Wow, that's really cool! Looks like this is a Haskell library: https://hackage.haskell.org/package/sbv-5.15/docs/Data-SBV.h...

Yes, it's really good! Makes it feasible to "sat first, ask questions later", especially for bit trickery.

Ah took me some time to understand why this works:

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.

Those are very convenient numbers for fizz buzz. Have you ever pulled out that trick in an interview?

Hmm, I wonder if this is more efficient than calculating mod 15 directly. It's basically something like

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

versus something like

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

Well, you're probably better of using the method in the article for the division, but this should work for small numbers.

Yes, just last week. I'm applying for a batch at Recurse Center as a sabbatical this fall, and wanted to put in something to show off a bit.