flat assembler
Message board for the users of flat assembler.
flat assembler
> Blog > Division by multiplication  revisited Goto page 1, 2 Next 
Author 

Back on the ASM Community message board, when there was a discussion about performing division by invariant divisor with a MUL instead of DIV, I looked to create an explanation of why this works through the theory of 2adic numbers. But recently, when I played with these ideas again, I realized that it is possible to analyze it just on the basis of pure integer arithmetic. I also made a new (at least to me) variant of the algorithm, which works correctly for every possible input and is easily provable.
This is how the algorithm looks implemented in fasm/fasmg scripting language: Code: ; Preparation stage:
BITNESS = 32 ; number of bits in word
DIVISOR = 3 ; a number to divide by
; convert DIVISOR into MAGIC and SHIFT constants:
SHIFT = bsr (DIVISOR1) + 1
MAGIC = ((1 shl SHIFT  DIVISOR) shl BITNESS) / DIVISOR + 1
; Division algorithm
NUMBER = 7FFFFFFFh ; an example number to divide
PRODUCT = NUMBER * MAGIC + NUMBER shl BITNESS
QUOTIENT = PRODUCT shr (SHIFT+BITNESS)
REMAINDER = ( (PRODUCT and (1 shl (SHIFT+BITNESS)  1)) * DIVISOR ) shr (SHIFT+BITNESS) The implementation in assembly language is a bit more complex. Read further below to learn why it has additional quirks compared to the classic ones: Code: ; Preparation stage:
mov rbx,[divisor]
mov rcx,rbx
dec rcx
bsr rcx,rcx
inc cl
mov rdx,1
shl rdx,cl
sub rdx,rbx
xor rax,rax
div rbx
inc rax
mov [magic],rax
dec cl
mov [shift_m1],cl Code: ; Division:
mov rbx,[number]
mov rax,[magic]
mov cl,[shift_m1]
mul rbx
add rdx,rbx
rcr rdx,1
shr rdx,cl
mov [quotient],rdx _____________________________________ And now the explanation. To analyze the problem purely from the point of view of the integer arithmetic, we can formulate it as follows. For a given invariant divisor D < 2^B (where B is the number of bits in word) we want to find a bit shift S and a "magic" number M such that: M*D = 2^(B+S) + R where 0 <= R <= 2^S If we succeed, then to get the quotient Q = N / D, we can do the following calculation: N*M = Q*2^(B+S) + L with 0 <= L < 2^(B+S) (this shows that we can drop the low B+S bits and take the high ones to get the Q) Let's take a look at the N*M*D product and substitute M*D from the definition of M: N*M*D = N*2^(B+S) + N*R Because N < 2^B and R <= 2^S, we have N*R < 2^(B+S), so the high bits of N*M*D contain original number N. If we look at the components of N*M, we can evaluate the same product in a different way: N*M*D = (Q*2^(B+S) + L)*D = Q*D*2^(B+S) + L*D now, as L < 2^(B+S), L*D < D*2^(B+S), therefore high bits of L*D contain a number smaller than D and we can represent it as: L*D = X*2^(B+S) + Y where X < D and Y < 2^(B+S) Now: N*M*D = Q*D*2^(B+S) + X*2^(B+S) + Y therefore high bits of N*M*D contain just Q*D + X and we already know that they also equal the original number N, so: N = Q*D + X and because X is positive and smaller than D, this proves that the quotient obtained this way is correct. Moreover, this demonstrates that high bits of L*D contain the remainder. _____________________________________ Now let's look at a classic example, this time for B=32: For D = 3 we can take M = 0AAAAAAABh and S = 1, then M*D = 200000001h = 2^(1+32) + 1 which fulfills the requirements, because 1 < 2 = 2^S. Therefore, to divide a number by 3, we can simply do this: Code: mov eax,[number]
mov edx,0AAAAAAABh
mul edx
shr edx,1
mov [quotient],edx Code: mov ecx,3
mov ebx,0
cmovc ebx,ecx
mul ecx
add edx,ebx
shr edx,1
mov [remainder],edx _____________________________________ For some divisors finding M and S such that R < 2^S may be troublesome. We can try the most straightforward route and just take smallest S such that D <= 2^S, then divide 2^(S+B) by D and increase it by one to cover the deficit caused by the truncation inherent in integer division: M = 2^(S+B) / D + 1 then 2^(S+B) < M*D <= 2^(S+B) + D, so the requirements are fulfilled, the only problem is that M is going to be larger than 2^B, so it no longer fits in a machine word. But, because we take S such that 2^S is the smallest power of two that is not smaller than D: 2^S / D = 1 and therefore M contains only a single bit set above the machine word. We can take M' = M  2^B and get the "magic" number that fits in the machine word. But when multiplying by it, we have to remember that we need to actually multiply by M' + 2^B, which is the same as multiplying by M' and then adding N to the high bits. This is how the algorithm I presented in the beginning works. The multiplier it computes is: (2^(S+B) / D + 1)  2^B = ((2^S  D) * 2^B) / D + 1 a then the product is adjusted by adding N*2^B before producing the quotient. This may produce carry, therefore the first of the right shifts that follows has to be RCL, to get the overflowing bit back into position. _____________________________________ If we prefer to stick to the classic algorithm with M fitting in the machine word, it may not always be possible to find the right combination of M and S as defined above. But we can also search for a pair such that: M*D = 2^(B+S)  R where R <= 2^S, as previously. What this modifies in the evaluations is that: N*M*D = N*2^(B+S)  N*R and since N*R again is smaller than 2^(B+S), the high bits of N*M*D this time are N1. So when we multiply by a number of this form, we actually compute (N1) / D instead of N / D. The classic algorithm deals with this by simply increasing the number by one before multiplication. For example, for D = 7 we can take M = 92492492h and S = 2. Then: M*D = 3FFFFFFFEh = 2^(2+32)  2 and 2 < 4 = 2^S. Therefore we can divide by 7 this way: Code: mov eax,[number]
mov edx,92492492h
inc eax
mul edx
shr edx,2
mov [quotient],edx 

06 Sep 2017, 21:31 

Tomasz Grysztar wrote: Therefore we can divide by 7 this way: Code: mov eax,[number]
mov edx,92492492h
inc eax
jz @f
mul edx
@@: shr edx,2
mov [quotient],edx The most common divisor 10 where M=CCC...CCD and S=3 can use the normal N*M >> S approach. 

07 Sep 2017, 11:28 

revolution wrote: If branches are cheap then we can add a single jz With this correction in mind, the only remaining advantage of the first algorithm is that it does not require any assumptions concerning divisors. But since the problematic divisors are the powers of two, we'd be better off handling them separately anyway. 

07 Sep 2017, 12:01 

Tomasz Grysztar wrote: It only gives the correct upper half of the product, though ... 

07 Sep 2017, 12:09 

revolution wrote:


07 Sep 2017, 12:21 

revolution wrote: Also if you know your numerator doesn't need it then leave out the jz. The requirement that R <= 2^S is only needed to show that N*R < 2^(S+B). So if we can get the same condition through the assumptions on N, it may enable new choices of S, even S = 0. So if we know that N is always small enough that N*D <= 2^B, we can take S = 0, M = 2^B / D + 1 (except when D is a power of 2, when we should drop "+ 1"). Then multiplying by M gives correct N/D in the high bits for all N < M. For example for D = 10, if we can safely assume that N is smaller than M = 2^32 / 10 + 1 = 1999999Ah, we can just multiply by it and take the high 32 bits as a result. 

07 Sep 2017, 12:21 

What about dividing numbers larger than the machine word? With DIV instruction it is possible to chain multiple divisions together, the remainder from the previous division is already in EDX/RDX ready to become the high word (smaller than the divisor, so overflow never happens) of the subsequent division in chain. With such chaining it is possible to divide arbitrarily large number by a small divisor.
The high word (EDX or RDX) must be smaller than the divisor, otherwise an overflow would occur (with DIV instruction it would be an exception 0). If we take V = bsr D, then: N < 2^(B+V+1), because high word of N is smaller than the divisor R < D < 2^(V+1) N*R < 2^(B+2*V+2) So we can take S = 2*V+2. Obviously M = 2^(B+S) / D + 1 is going to be larger than machine word, so the the multiplication is no longer going to be simple. For D = 10, we get S = 8 and M = 2^40 / 10 + 1 = 199999999Ah. This led me to the following 32bit routine that can replace "DIV 10" (without checking for overflow, input EDX must be smaller than 10): Code: div10:
; this trashes EBX and ECX, you may need to preserve them
push eax edx
imul edx,19h
mov ecx,edx
mov edx,09999999Ah
mov ebx,eax
mul edx
mov eax,ebx
mov ebx,edx
mov edx,19h
mul edx
add ebx,eax
adc ecx,edx
pop eax
mov edx,09999999Ah
mul edx
add eax,ebx
adc ecx,edx
shrd eax,ecx,8
pop edx
imul ecx,eax,10
sub edx,ecx
retn 

08 Sep 2017, 15:05 

Since this is a thread on division, I though I would collect my thoughts on multiprecision division. Throughout this post, let be the word size of the machine  so usually , but any integer will do. Given wordsized positive integers , the x86 architecture has an instruction div for calculating
An exception is raised if this quotient is , i.e. not word sized. The goal of division is to calculate for multiword integers and . When is wordsized, this can be calculated by stitching together div's, so the problem is to solve the case when is bigger than a word. In the grade school algorithm, the words of the quotient are produced from most significant to least significant. This means that the grade school algorithm is only dependent upon a solution to the following problem: Given big integers and with and , calculate the wordsized integer . Now, by shifting and , this problem can also be represented as calculating where and are integers with and , and and are real numbers with . Fact: With the above assumptions, it holds that . This means that the difficult (n+1)byn word division can be reduced to the simpler 3by2 division . If the latter does not match the former, is must be over by only one. In the grade school algorithm, this corresponds to overestimating a quotient term, and thus having to 'fix it up'. The point is that if the quotient term is estimate incorrectly, it is overestimated and only overestimated by one. Proof: Write where and . Then, It is easy to see that and some manipulation of the assumption shows that This last quantity is greater than by the assumptions on . So now the grade school division has been reduced to the division problem with the assumptions that the and are wordsized and that . This will be covered in a future post. 

09 Sep 2017, 09:48 

In the previous post, large integer division was reduced to a sequence of 3by2 divisions followed by nby1 multiplysubtract operations, in true grade school fashion. The quotient of the 3by2 division was between and inclusive. Since the quotient terms for the grade school division are by definition wordsized, we need to calculate
where it is assumed that this floor is and of course all of the previous assumptions, including . If you follow the proof given in the previous post for this problem you will see that However, care is required as the division can overflow a word. I propose that the following algorithm computes . Code: if B2 < A1
(q, r) = div(B2:B1, A1)
(lo, hi) = mul(q, A0)
T0 = sub(B0, lo)
T1 = sbb(r, hi)
T2 = sbb(0, 0) ; copy carry flag across T2
if T2 <> 0
q = dec(q)
T0 = add(T0, A0)
T1 = adc(T1, A1)
T2 = adc(T2, 0)
if T2 <> 0
q = dec(q)
T0 = add(T0, A0) ; only calculated
T1 = adc(T1, A1) ; for the asserts
T2 = adc(T2, 0) ; below
end if
end if
assert( T2 == 0) ; T2:T1:T0 is supposed
assert( T1 * X + T0 < A2 * X + A1) ; to contain the remainder
return q
else
q = X  1
T0 = B0
T1 = sub(B1, A0)
T2 = sbb(B2, A1)
T3 = sbb(0, 0)
T0 = add(T0, A0)
T1 = adc(T1, A1)
T2 = adc(T2, 0)
T3 = adc(T3, 0)
if T3 <> 0
q = dec(q)
end if
return q
end if I do not know at this time why at most one fixup is required in the "else" block. Of course two fixups are required in general by the above inequality. Also, notice that this algorithm contains a call to div with a fixed divisor for each quotient word produced in the whole grade school algorithm. Thomasz would probably want to change this invariant division into some multiplications. 

09 Sep 2017, 11:07 

When I was doing this stuff a few years ago I found that Barrett and Montgomery were the most used two methods for division and modulo for multiprecision operations. There are also whole books dedicated to this subject, and many PhD theses that attempt to improve the efficiency and performance in specific CPUs.


09 Sep 2017, 11:59 

tthsqe wrote: Since this is a thread on division, I though I would collect my thoughts on multiprecision division. Also, what I wanted to focus on here is the analysis of these tricks purely on the grounds of integer arithmetic, with no event slightest mention of real (or even padic) numbers. This is my personal hobby. 

09 Sep 2017, 13:28 

I don't think it needs to be split. I am just giving you more to think about.
What you have been doing is the division of singe words (like the c "/" operator). You can interpret what I am saying as request for a division of two words by one word (like the x86 "div" instruction) by using multiplications with precomputed numbers. 

09 Sep 2017, 14:50 

Back to the beginning: the idea behind the division by multiplication can be presented as follows.
If we want to divide integer N by divisor D, what we can try is to find number such that and: , . The same product can also be written as: and because : , so: , , . Therefore and thus which gives us correct division with remainder, since . We get the result Q simply by taking the high bits of X. We can also get the remainder R by taking high bits of LD. Now, to find the number X we can multiply N by socalled magic multiplier M, which is a number such that: , and then: which gives the correct result as long as . All the algorithms presented above made this work by finding n large enough that this inequality would be fulfilled for all N in the required range. But there is still another approach we could try. When we compute M, we also know P, and we can use this knowledge to try to adjust the result even if NP ends up being large enough to meddle with the high bits of NMD. We can compute NP and split into high and low part: and then try to adjust for the error by taking: . This leads to: Now, keeping in mind that: we also need to ensure that: . Then the worst that may happen is that the high bits of N'MD might contain N  1 instead of N. This can be easily fixed after computing the remainder  if remainder ends up being larger than D we can subtract D from it and increase the quotient by one. Let's try it on an example, with D = 10, n = 32. Then M=1999999Ah and P = 4: To simulate a DIV instruction we need to operate on range: With this in mind we have: so and finally: This is far from being dangerous for the high bits of N'MD. It is time for an implementation: Code: div10:
; cmp edx,10
; jae int0
push ebx ecx
push eax
mov ebx,eax
mov ecx,edx
shld edx,eax,2
sub ebx,edx
sbb ecx,0
mov eax,ebx
mov ebx,1999999Ah
mul ebx
mov eax,ecx
imul eax,ebx
add eax,edx
pop edx
imul ecx,eax,10
sub edx,ecx
cmp edx,10
jb .done
sub edx,10
inc eax
.done:
pop ecx ebx
retn For 64bit version just replace all 32bit registers with their 64bit counterparts and the constant with 199999999999999Ah (P is the same). 

09 Sep 2017, 14:53 

tthsqe wrote: You can interpret what I am saying as request for a division of two words by one word (like the x86 "div" instruction) by using multiplications with precomputed numbers. . But when dividing large integer by another large integer, you usually need to divide by a normalized divisor, such that: . I have not yet tried to solve this variant with my integerbased approach. When a divisor is normalized, the algorithm I presented first has the shift equal to the size of machine word (the assembly implementation I provided does not work then, because SHRD sees such shift as zero), so the result can be taken straight from the third word after multiplication. Also, multiplying the lowest word seems pointless, because it can affect the result by no more than one, and this can be corrected when computing remainder, as in my previous post. However to extend the range of algorithm, some kind of compensation for the NP term would be needed. I don't know if I can solve this variant with integer arithmetic as nicely as the previous ones. 

09 Sep 2017, 14:58 

Tomasz Grysztar wrote: This bring us to the following implementation of the classic algorithm: I'm trying to write a compiletime tool to help with generating magic numbers and shifts in the most optimal way possible at compiletime (so I don't care about the runtime divisors), but they have to be generic (the +1 condition is an extra insn other than mul/shift). For now I use this algorithm for fullwidth multiplies, and your first algorithm for smaller muls, because then the mask with width+1 bits will fit and I won't have to do the +1 (inc) in some cases (just a mul and shift). But as you said: Tomasz Grysztar wrote: Note that both the algorithms given use the large values of S even if there may exist a smaller S for which the requirements would be fulfilled. An example for B=32 is number 641, where the algorithm works for S = 0 and M = 663D81h: Also, is there a better way to implement such divisions by constants when the numbers are smaller than fullwidth? (e.g. 16bit or 8bit numbers involved, or even 32bit if x64). I mean, I can obviously use the entire register for an 16bit operation so was wondering if there's a better way (I currently use your first method with the mask overflowing into 16th bit for example, for 16bit number, and do a 32bit mul). What I consider optimal: least amount of extra instructions other than mul. (and yes, I will also implement this tool in C++ as a compiletime metaprogramming thing not just an asm tool, I just found out GCC can be quite retarded sometimes, so...) 

28 Feb 2018, 15:56 

Furs wrote: Hi, I'm wondering, does this algorithm fail for anything other than power of 2s? (those can be easily handled) Furs wrote: Is there a way to reliably handle such cases? It doesn't really matter if it's a bit "expensive" to detect (loops?) since it is a tool done at compiletime, not runtime, which I don't care about. i.e. I just want it to spill out a magic, a shift (which can be 0 if it's more optimal) and an add (which is 0 most of the time). Furs wrote: Also, is there a better way to implement such divisions by constants when the numbers are smaller than fullwidth? (e.g. 16bit or 8bit numbers involved, or even 32bit if x64). I mean, I can obviously use the entire register for an 16bit operation so was wondering if there's a better way (I currently use your first method with the mask overflowing into 16th bit for example, for 16bit number, and do a 32bit mul). Hence to divide 16bit N by 16bit D, you can use M = 2^32 div D, do a 32bit multiplication and take DX as a result (high word of EDX is always going to be zero). Last edited by Tomasz Grysztar on 01 Mar 2018, 10:54; edited 1 time in total 

28 Feb 2018, 16:54 

I'm a bit confused. So I need to iterate all shift values, but how can I check if the condition holds if I don't know the value of N? Should I assume the highest value possible for it that it can have? I mean, N is a runtime value, only the divisor is at compile time.
And also if I divide 2^32 by 641, the remainder is 640, not 1, so how come it works? Tomasz Grysztar wrote: What limits the range where the algorithm works is the already mentioned constraint. If you choose to not use shift and take n equal to the size of register, you can safely work with values of half that size. For example if n = 32 and both N and D are 16bit numbers, then also P fits into 16 bits (since it is a remainder, smaller than D) and the inequality holds. I think I understand it now. Sorry for being slow, I'm not that good at mathspeak. EDIT: I'm guessing P is not the remainder, but rather Dremainder, am I right? I looked before at your example with 10, and you said P = 4, which fits just as well as on 641. (i.e. 106 = 4, and 641640 = 1) 

28 Feb 2018, 18:16 

I implement this from
http://www.flounder.com/multiplicative_inverse.htm I test with some values and work well, of course for power of two I use the shift transformation. I'm implementing now transform multiplication with shift and adds, but search the optimal combination is hard, few web pages with information, I try to dig in source code of gcc but not found this. 

28 Feb 2018, 21:19 

pabloreda wrote: I'm implementing now transform multiplication with shift and adds, but search the optimal combination is hard, few web pages with information, I try to dig in source code of gcc but not found this. Anyway, this is kind of what I have right now, in pseudocode (I have it in C++ and inline asm (for rcr) but I figure it's not so easy to understand): Code: ; 'max' is the max value of the numerator N (the output depends on the programmer knowing this in advance)
; N is the numerator
; D is the divisor
if(D is power of 2) return N / D; // compiler handles this here
for(shift in 0 ... 31) {
temp = 2^32 << shift;
if(max*(D  temp mod D) < temp) break;
}
if(shift < 32 && (magic, which is (2^32 << shift) / D, fits in 32 bits)
return (N*(magic+1) >> 32) >> shift;
// check if we can't use 'inc' method
if(max == 0xFFFFFFFF)
{
// use Tomasz's method
shift = bsr(D1) + 1;
magic = (((2^32)b) << 32)/b + 1;
return (((N*magic >> 32) + N) rcr 1) >> (shift1);
}
else
{
shift = bsr(D);
magic = (2^32 << shift) / b;
return (N+1)*magic >> 32 >> shift;
} 

28 Feb 2018, 21:37 

Goto page 1, 2 Next < Last Thread  Next Thread > 
Forum Rules:

Copyright © 20042018, Tomasz Grysztar.
Powered by rwasa.