How can I reverse optimized integer division/modulo by constant operations?

Question

When compiling a division or modulo by a constant, my compiler (LLVM GCC) generates a series of instructions that I don't understand.

When I compile the following minimal examples:

int mod7(int x) {
    return x % 7;
}

int div7(int x) {
    return x / 7;
}

The following code is generated:

_mod7:
    push   rbp
    mov    rbp,rsp

    mov    ecx,0x92492493
    mov    eax,edi
    imul   ecx
    add    edx,edi
    mov    eax,edx
    shr    eax,0x1f
    sar    edx,0x2
    add    edx,eax
    imul   ecx,edx,0x7
    mov    eax,edi
    sub    eax,ecx

    pop    rbp
    ret    


_div7:
    push   rbp
    mov    rbp,rsp

    mov    ecx,0x92492493
    mov    eax,edi
    imul   ecx
    add    edx,edi
    mov    ecx,edx
    shr    ecx,0x1f
    sar    edx,0x2
    mov    eax,edx
    add    eax,ecx

    pop    rbp
    ret

• How is this mathematically equivalent, and where do the constants come from?
• What's the easiest way to turn the assembly back in to C (for arbitrary constants on the right-hand side)?
• How could a tool, such as a decompiler or analysing disassembler, automate this process?

Answer 1

First

Unfortunately we don't seem to have MathJax turned on in this stackexchange so the math parts below are pretty horribly formatted. I'm also far from a mathematician so the notation may be off in some places.

Understanding the magic number and code

The goal of the code above is to rewrite a division into a multiplication because division takes more clock cycles than a multiplication. It's in the area of about twice as many cycles, depending very much on CPU. So we need to find a nice branchless way of doing that. If we branch we're very likely to lose to simply doing division.

One way is to simply realize that division is the same as multiplication with the inverse of the number, i.e . The problem is that is a pretty poor number to store as an integer. So we need to multiply both the divisor and dividend by some number. Since we're operating on 32-bit numbers and we get multiplication results in 64-bit numbers we get the best precision with and we also avoid overflow issues. So we basically get . Now that fractional part is what causes us issues because it will cause rounding errors.

So let's try to formalize this:

Where is our multiplicand, e.g , or really any number but works very well with our register sizes as we can simply discard the lower 32-bit register. is the number you must add to make evenly divisible by . is the number we wish to divide.

We can rewrite the equation above, as

Which illustrates the point that we have our dividend divided by our divisor and then an error factor of .

Studying our original equation of it's clear that we can affect very little. needs to be a power of 2, can't be too large or we risk an overflow and can't be too small as it has a direct negative effect on our error factor . directly depends on and .

So let's try which gives a maximum error fraction of with the maximum value of being , so , unfortunately this is not less than so we can get rounding errors.

We'll increase the exponent of to , which gives , maximum error fraction which is less than . This means that our multiplicand is which is not less than or equal to the maximum signed value we can store in a 32-bit register (). So we instead make the multiplicand . As a side note, thanks to the magic of two's complement when we subtract the number is which is when interpreted as an unsigned number. But we're doing signed arithmetic here. So we need to fix the final expression by adding . This also only solves the problem for , for negative numbers we will be off by 1 so we need to add 1 if we have a negative number.

That's the explanation for the constant in the multiplication and how to arrive at it. Now let's look at the code:

; Load -1840700269
mov    ecx,0x92492493

; Load n
mov    eax,edi

; n * -1840700269
imul   ecx

; add n to compensate for 2^32 subtraction
add    edx,edi

; check the sign bit of our result
mov    ecx,edx
shr    ecx,0x1f

; divide by 2^2 to compensate for us using y=2^34 instead of 2^32
sar    edx,0x2

mov    eax,edx
; add the value of the sign bit to the final result
add    eax,ecx

Calculating divisor from magic number and code

I have not proven this mathematically, however if you want to recover the divisor from an assembly dump such as the one you showed we can do some simple mental excercises. First we need to realize that the following holds

Where is the adjustment we made in order to bring the value into the range of a 32-bit value. From the code we can devise the following, the right shift by two means that we have , , , is unknown. This means that we're missing one variable in order to perform a perfect solution. However the effect of if negligible as its purpose is to bring the divisor as close to its integer value as possible. This means that the solution can be found by solving

Another example with divisor 31337 which has the multiplicand magic number 140346763 and right shifts 10 bits.

Finally

For a complete mathematical breakdown of how this works, including all the appropriate proofs and algorithms for calculating the magic numbers, shifts and adds, see Hacker's Delight, chapter 10-3.

Answer 2

Here's a late response. The Reko decompiler recovers the integer divisors by performing a divide and conquer search using mediants.

Reko starts by recognizing the pattern where the high word of a 64-bit product (r * c) is used. The constant multiplier c is divided by 2^32 to yield a double precision floating point number between 0.0 and 1.0. Starting with the rational numbers 0 / 0 and 1 / 1, Reko computes a sequence of mediants that brackets the floating point number. From this sequence of mediants it chooses the rational number that comes closest to the floating point number and returns it.

The code is not fully tested yet -- I haven't had a chance to work with negative numbers yet, for one, but seems to give reasonable results. The code is here if you're curious: https://github.com/uxmal/reko/blob/master/src/Decompiler/Analysis/ConstDivisionImplementedByMultiplication.cs

Answer 3

This paper might be of interest: Division by invariant multiplication.

Bumped into this here.