glen said:
Terje Mathisen wrote:
(snip)
That is nice. I probably would have guessed that there would be some
cases where it rounded wrong, but that would have been a guess. I
would have guessed that it would take at least one more bit than the
dividend length to get it right.
That's correct: You do need a 33-bit reciprocal, rounded up.
This means that for that half of divisors where the reciprocal's 33rd
bit is 1, you get to round it up to a 32-bit reciprocal (with trailing
zeroes), but for the other half you need some more work:
Either you must add half the input value once more to the low half of
the product, then propagate any carries, or you can take advantage of
the fact that the top reciprocal bit is always 1, so you can multiply by
the next 32 bits, then add the missing part afterwards. Both of these
suffer from the need to temporarily work with more than 32 bits:
First method
mov eax,[reciprocal_rounded_down] ; First 32 bits of reciprocal
mul ebx ; Value to be divided
shr ebx
add eax,ebx
adc edx,0
mov cl,[shift]
shr edx,cl ; EDX has exact result
Second method:
mov eax,[reciprocal_without_top_bit]
mul ebx
add edx,ebx ; Can overflow into 33 bits!
rcr edx,1
mov cl,[shift]
shr edx,cl ; EDX has exact result
Agner Fog & I corresponded on this algorithm for a while, he discovered
a much better approach to use in the case of too-long reciprocals:
Instead of faking a 33-bit reciprocal, he determined that it is OK to
adjust the input value instead, and use a rounded-down 32-bit reciprocal:
mov eax,[reciprocal_rounded_down]
add ebx,1 ; Value to be divided, adjusted up. (Can overflow!!!)
jc skip_mul
mul ebx
skip_mul:
mov cl,[shift]
shr edx,cl ; EDX has exact result
If the cost of that branch can be high, either because maximal inputs
are relatively common, or because your code is running out of branch
predictor entries, then you can make it branchless like this:
mov eax,[reciprocal_rounded_down]
inc ebx ; Value to be divided, adjusted up. (Can overflow!!!)
mul ebx
test ebx,ebx
cmovz edx,[reciprocal_rounded_down]
mov cl,[shift]
shr edx,cl ; EDX has exact result
The best choice for a semi-constant division would be an algorithm that
always works, right? It turns out that you can do this in various ways,
but the trailing 32-bit reciprocal + extra add is probably the easiest:
Works for any divisor, calculate a 33-bit accurate reciprocal, where the
33rd bit is rounded up, then remove the top bit:
mov eax,[reciprocal_without_top_bit]
mul ebx
add edx,ebx ; Can overflow into 33 bits!
rcr edx,1
mov cl,[shift]
shr edx,cl ; EDX has exact result
Now, will it still work for 64 bits?
Yes, as long as you follow the same rules!
In fact, on a 64-bit system you can use the same approach even in C as
long as you have some kind of compiler intrinsic to give you access to
the high half of a 64x64->128 bit multiplication. You can do it even if
you don't have a carry flag to temporarily give you an extra bit.
I.e. this works on Alpha and most others.
uint64_t reciprocal_without_top_bit;
unsigned shift;
void setup_divisor(uint64_t divisor);
inline
uint64_t semi_fixed_divide(uint64_t n)
{
uint64_t prod = MULH(n, reciprocal_without_top_bit);
prod += n; // Adjust for missing top bit
int carry = prod < n; // Overflow? 0 or 1
prod = (prod >> 1) | (carry << 63);
prod >>= shift;
return prod;
}
Terje
