Division Algorithm for microcontrollers

(c) 1999, 2003 EME Systems, Berkeley CA U.S.A.
Tracy Allen
<stamp index> <home>

Contents (updated 11/19/06)

 

Algorithm for multi-byte division on microcontrollers:

This note discusses the process of division on microcontrollers, not specifically the BASIC Stamp. Assuming we are dividing one positive integer N by another D, both possibly multiple bytes. We will solve for quotient Q and remainder R in

        N = Q*D+R, where 0<=R<D, Q an integer

or to write it another way:

        N/D = Q  +  R/D   where 0<=R/D<1.

Steps 1 to 4 solve this equation to find the integer part. Steps 5 to 8 extend the result to expand the remainder R/D as a fixed point fraction.

This is followed by a discussion of the theory, and the interpretation of fixed point fractions in binary and decimal forms.

 

N  the integer numerator, multiple bytes possible
kw the register size in bits (eg 64 bits for 8 bytes)
kn the number of significant bits in the numerator, to be determined in real time.
D the integer denominator, multiple bytes
Q the integer quotient, multiple bytes (could be as large as N if D=1)
R the integer remainder, multiple bytes (allow width up to D*2 for computation)
X a working register, same size as R
i a loop index, byte

Assume we have the following subroutines among other minor ones:

X=:R-D multibyte subtraction ends with carry set if D>R.
N=:N<<1 or D=:D<<1 or Q=:Q<<1 multibyte rotate through carry

Steps in pseudocode:

1) if N=0, exit to special case, numerator=0
if D=0, exit to special case, denominator=0
account for negative values, if needed. This routine works only on + values
if N<D, set Q=0, R=N, and go to step 4, there is no integer part.
2) determine kn, the number of significant bits in N, e.g. 10111: kn=5
N:=N<<(kw-kn) left justify the numerator in the kw bit field
R:=0 initial remainder
Q:=0 initial quotient
i:=kn+1 initial loop index
3) i:=i-1
if i=0 then go to step 4
N:=N<<1 shift one bit left, lsb:=0, carry:=msb
R:=R<<1 shift carry into lsb of R.
X:=R-D subtract sets carry if D>R
complement carry carry is set if R>=D
Q:=Q<<1 shift carry from the subtraction into lsb of Q
go back to 3 if carry is set from the subtraction step
R=X replace R with R-D
goto 3
4) Now have integer quotient Q and remainder R. Stop if that is all we need.
We have, N = Q*D+R, 0<=R<D
or written another way
N/D = Q + R/D 0<=R/D<1
5) now to find the binary representation of the fractional part, R/D, to be carried out in powers of 1/2^m to z steps,.
F:=0 initial quotient (in approximation R/D=F/(2^z)
S:=R initial remainder (from above)
D denominator (from above)
i:=z+1 initial index
6) i:=i-1
if i=0 then goto step (7)
S:=S<<1 shift S left one bit, 0 into lsb
X:=S-D subtract sets carry if D>S
complement carry carry is set if S>=D
F:=F<<1 shift carry into lsb of F; least sig bit set if S>=D
go back to 6 if carry is set
from subtraction step
  S=X replace S with (S-D)
goto 6
7) Now have fraction F, as approximation to R/D=F/(2^z)
You can look at this as the z binary digits to the right of the binary "decimal" point.

8) convert F to decimal fraction, if desired.
The fraction we have approximated is R/D, and we have done that by
finding an approximation that has a denominator with a power of two.
Ususally the power of two will be some multiple of the word size,
but it can be any power of two that is convenient, by concluding
the above algorithm at the desired number of bits.
Now we may want to convert that to a decimal fraction, with a power
of 10 in the denominator:

R/D = F/(2^z) = F'/(10^y)
F':= F*(10^y)/(2^z)
original value of F is multiplied times the binaray representation
of 10^y, and then the result is divided by 2^z (by shifting or by
taking the upper bits of the double product. This is like the **
operator in the BASIC Stamp.

Step 3 can be shortened at the beginning by accounting for leading zeros in the quotient, i.e., the remainder must have at least the same number of bits as the denominator, before a 1 will first appear in the quotient.

Note that step 6 is almost the same as step 3.   The only difference is that there is no longer an integer part, N in step 3.  By the end of step 3, the integer part has been reduced to zero.  At the end of step 6, there still may be a remainder in S.   The algorithm could be continued to kick out additional binary digits to the right of the binary "decimal" point.   The fraction terminates eventually only if the fraction can be written as exact powers of (1/2),   The denominator is a power of two only.


Theory

The algorithm in step 3 generates the polynomial:

   Qm = [f(m)*2^(m) + f(m-1)*2^(m-1) + ... + f(1)*2 + f(0)

The most significant bits of the numerator are shifted into R one at a time, until the value of R is large enough to be divided by D. At that point a 1 is shifted into the quotient, and D is subtracted from R and a new bit of the numerator is shifted in at the right. The process continues until no more bits are left in the numerator. If D is greater than the remainder R at a step in the process, a zero is shifted into the quotient.

If a remainder is left at the end of the division, the process can continue (step 5 to 8) to generate a fixed point fraction.

If the initial values of N and D are already contain a decimal point, some modification is necessary. Of course, it is not really a "decimal" point, but a "binary" point. See below for discussion of the meaning of fixed point fractions as represented in the binary or decimal systems. The point in D should be shifted right to the right far enough to clear the fraction, and the point in N must be shifted a like amount. If there is still a fractional part in N, a further shift can be carried out to get rid of that too, so that the division can take place using large integers only. Or the fixed point fractional part of the numerator can be rolled into the calculation of step 6 above with only slightly more complication, like step 3.

Here is more explanation of the fractional part. Confusion arises about the meaning of fractions on a binary computer, compared to decimal fractions.

We have the initial computation:

   N = Q * D + R

The quotient Q is the integer part of N/D, and the remainder also is an integer R, that satisfies the inequality

   0<=R<D.

The equation can be written:

  N/D = Q + R/D    with 0<=R/D<1

The fractional remainder, R/D is a number less than one. Let us continue to expand R/D as a polynomial. Again this is easiest to implement on a binary computer by using the binary arithmetic.

The numerator in the approximation, Fm, is a polynomial in inverse powers of two, with binary coefficients, that are "generated" by the algorithm. The binary fraction is:

   R/D = f(0)*1/2 + f(1)*1/4 + ... + f(m-1)*1/2^(m-1) + f(m)/(2^m)

and pulling out the factor 1/2^m from the denominator, we have

   R/D = [f(0)*2^(m-1) + f(1)*2^(m-2) + ... + f(m-1)*2 + f(m)]/2^m

= Fn/2^m

The algorithm above generates the polynomial Fn.


Interpretation of the fixed point fraction

A word about interpretation of the fixed point fraction: Given numerator over denominator, R/D, and m-step binary algorithm, we are really forming the approximation

  R/D = Fn/(2^m)   (remainder truncated)

That is, we are approximating N/D by a fraction where the new denominator is a power of 2, e.g.

   4567/6789 =~ 2889249615/4294967296   to this precision

=~ 44086/65536 to earlier step of precision

If we use a decimal algorithm it is the same idea, but the denominator is a familiar power of 10, e.g.,

   4567/6789 =~ 672705/1000000    to 6 places

=~ 672/1000 to 3 places, earlier step of approximation

The first point is, that if we use the binary algorithm, which is very fast and easy to implement on a binary computer, we will find a result that is not a decimal (in the sense that a decimal fraction has a power of 10 in the denominator, not a power of 2). To convert the result of the binary algorithm to decimal form, we have in effect to solve the equation:

    y/10^k = x/2^m

given x, find y, which simply requires multiplying the result of the binary division times the desired decimal denominator:

    y = x * (10^k) /(2^m)

This is a simple multiplication, followed by a shift, or a "take the upper bytes" trick like the ** operator on the BASIC Stamp. If this is an intermediate calculation in a longer procedure, it may be possible to postpone the renormalization until the end.

Renormalization is not a consideration while calculating the integer part of a quotient. The quotient Q in N/D = Q + R/d is a simple integer, not a ratio.

Now, the algorithm for determining the fixed point fraction. The following may confuse, or clarify. Just an talk-it-through exercise.

At each step in the division, the result is a quotient and a remainder.

         N  * (2^m) =   Fm  *   D  + Rm

4567 * 65536 = 44086 * 6789 + 3058 <--- example

Note that if you stop at that step and truncate the remainder, you get the approximatione:

    N/D = Qn/65536  = 44086/65536 (remainder truncated)

The next step is multiply both sides by 2 (shift left), and divide again by D. Let n=m+1 be this next subscript:

          N * (2^n) /D  =   Qm*2    + Rm*2/D

All the powers in the polynomial Qm are multiplied by 2, that is, each exponent increases by 1. The quotient Rm*2/D has an integer part, which is either 0 or 1, and a fractional part, Rn/D. The integer part of Rm*2/D (either q0=0 or q0=1) is now the coefficient of 2^0 in the new polynomial, Qn. The new remainder is either Rn=Rm (if q0=0) or Rn=Rm-D (if q0=1). The new Rn will be in the interval D>Rn>=0.

          N *  2^n =    Qn   *  D   + Rn
e.g. 4567* 131072 = 88172 * 6789 + 6116 <--- example

You could stop there and have the better approximation:

     4567/6789 = 88172/131072

Or you could perform exactly the same step again and again to improve the size, and precision, of the implied denominator.


<top> <index> <home> logo < mailto:info@emesystems.com >