janmr blog

Basic Multiple-Precision Long Division

We consider the task of dividing a positive integer uu by another positive integer vv, thus obtaining a quotient q=u/vq=\lfloor u/v \rfloor and a remainder rr such that u=qv+ru = q v + r with 0r<v0 \leq r < v.

The method presented here is based on The Classical Algorithms, Section 4.3.1, of The Art of Computer Programming, Volume 2, by Donald E. Knuth. The material is quite theory-heavy and if you are just looking for the main algorithm, you can skip to the bottom and Algorithm L.

We represent the numbers using radix b2b \geq 2 and set

u=(um1u1u0)bandv=(vn1v1v0)b  ,u = (u_{m-1} \ldots u_1 u_0)_b \quad \text{and} \quad v = (v_{n-1} \ldots v_1 v_0)_b \; ,

so uu is an mm-digit number and vv is an nn-digit number (see previous post for more details on representing multiple-precision numbers).

Two special cases are easily dealt with:

  • If m<nm < n then u<vu < v and so q=0q = 0 and r=ur = u is the simple answer.
  • If n=1n = 1 then vv is just a single digit and we use a short division algorithm instead.

So in the following we assume that mn>1m \geq n > 1.

We will approach the division algorithm from a top-level point of view. It is actually just a formalization of the well-known pencil-and-paper method:

Algorithm G. Given u=(um1u1u0)bu = (u_{m-1} \ldots u_1 u_0)_b, um10u_{m-1} \neq 0 and v=(vn1v1v0)bv = (v_{n-1} \ldots v_1 v_0)_b, vn10v_{n-1} \neq 0, with mn>0m \geq n > 0, this algorithm outlines how to compute the quotient q=(qmnq1q0)b=u/vq = (q_{m-n} \ldots q_1 q_0)_b = \lfloor u/v \rfloor (we may have qmn=0q_{m-n} = 0 in which case qmn10q_{m-n-1} \neq 0 if m>nm > n) and the remainder rr such that u=qv+ru = q v + r, 0r<v0 \leq r < v.

  • G1. u(mn+1)(0um1u1u0)bu^{(m-n+1)} \leftarrow (0 u_{m-1} \ldots u_1 u_0)_b.
  • G2. kmnk \leftarrow m-n.
  • G3. qk(uk+n(k+1)uk+1(k+1)uk(k+1))b/vq_k \leftarrow \left\lfloor (u^{(k+1)}_{k+n} \ldots u^{(k+1)}_{k+1} u^{(k+1)}_k)_b / v \right\rfloor.
  • G4. Set u(k)u(k+1)qkbkvu^{(k)} \leftarrow u^{(k+1)} - q_k b^k v or, equivalently,

(uk+n(k)uk+1(k)uk(k))b(uk+n(k+1)uk+1(k+1)uk(k+1))bqkv,(uk1(k)u1(k)u0(k))b(uk1(k+1)u1(k+1)u0(k+1))b.\begin{aligned} (u^{(k)}_{k+n} \ldots u^{(k)}_{k+1} u^{(k)}_k)_b &\leftarrow (u^{(k+1)}_{k+n} \ldots u^{(k+1)}_{k+1} u^{(k+1)}_k)_b - q_k v, \\ (u^{(k)}_{k-1} \ldots u^{(k)}_1 u^{(k)}_0)_b &\leftarrow (u^{(k+1)}_{k-1} \ldots u^{(k+1)}_1 u^{(k+1)}_0)_b. \end{aligned}

  • G5. If k=0k=0 then set ru(0)r \leftarrow u^{(0)} and exit. Otherwise set kk1k \leftarrow k-1 and go to step G3.

An essential invariant of this algorithm is

(uk+n1(k)uk+1(k)uk(k))b<vfork=0,1,,mn+1.(u^{(k)}_{k+n-1} \ldots u^{(k)}_{k+1} u^{(k)}_k)_b < v \quad \text{for} \quad k=0, 1, \ldots, m-n+1.

This can be seen as follows. For k=mn+1k=m-n+1 the invariant is ensured by introducing a zero as the most significant digit of u(mn+1)u^{(m-n+1)} in step G1. For k=0,1,,mnk=0,1,\ldots,m-n we see from steps G3 and G4 that (uk+n(k)uk+1(k)uk(k))b=(uk+n(k+1)uk+1(k+1)uk(k+1))b mod v(u^{(k)}_{k+n} \ldots u^{(k)}_{k+1} u^{(k)}_k)_b = (u^{(k+1)}_{k+n} \ldots u^{(k+1)}_{k+1} u^{(k+1)}_k)_b \text{ mod } v and the inequality follows.

Note that the invariant implies that uk+n(k)=0u^{(k)}_{k+n}=0 for k=0,1,,mnk=0, 1, \ldots, m-n. Furthermore we have that

(uk+n(k+1)uk+1(k+1)uk(k+1))b=(uk+n(k+1)uk+1(k+1))bb+uk(k+1)(v1)b+(b1)=vb1(u^{(k+1)}_{k+n} \ldots u^{(k+1)}_{k+1} u^{(k+1)}_k)_b = (u^{(k+1)}_{k+n} \ldots u^{(k+1)}_{k+1})_b \cdot b + u^{(k+1)}_k \leq (v-1) b + (b-1) = v b - 1

from which we see that the quotients qkq_k computed in step G3 are non-negative and smaller than bb, as they should be.

Finally, we can verify that the algorithm computes what we intended. We have

r=u(0)=u(1)q0b0v=u(2)q1b1vq0b0v==u(mn+1)(qmnbmn++q0b0)v=uqv.\begin{aligned} r &= u^{(0)} = u^{(1)} - q_0 b^0 v = u^{(2)} - q_1 b^1 v - q_0 b^0 v = \ldots \\ &= u^{(m-n+1)} - (q_{m-n} b^{m-n} + \cdots + q_0 b^0) v = u - q v. \end{aligned}

Now for some practical aspects. Note first that all of the u(k)u^{(k)} variables can in fact be represented by a single variable and simply overwrite its digits along the way – thus ending up with the remainder. Note also that any of the remainder's digits may be zero.

Finally, how do we compute the quotient in step G3? That is in fact the central part of the division algorithm and is the subject of the rest of this post.

Let us now consider computing the quotient in step G3 in the case n>1n > 1. We therefore assume u=(unu1u0)bu = (u_n \ldots u_1 u_0)_b, u<bn+1u < b^{n+1}, and v=(vn1v1v0)bv = (v_{n-1} \ldots v_1 v_0)_b, bn1v<bnb^{n-1} \leq v < b^n, with n2n \geq 2 and 0u/v<b0 \leq \lfloor u/v \rfloor < b.

We wish to compute q=u/vq = \lfloor u/v \rfloor as fast as possible. How good is a 'first order' approximation, where we use just the two most-significant digits of uu and the most-significant digit of vv: (unb+un1)/vn1(u_n b + u_{n-1})/v_{n-1}? First of all, if un=vn1u_n = v_{n-1} this quantity equals bb and we know that qb1q \leq b-1 by assumption, so let us therefore study

q^=min(unb+un1vn1,b1)\hat{q} = \text{min} \left( \left\lfloor \frac{u_n b + u_{n-1}}{v_{n-1}} \right\rfloor, b-1 \right)

This approximate quotient is never too small, as the following theorem states.

Theorem 1. With q^\hat{q} as defined above we have qq^q \leq \hat{q}.

Proof. If q^=b1\hat{q}=b-1 then since qb1q \leq b-1 by assumption, the statement is true.

Assume then that q^=(unb+un1)/vn1\hat{q} = \lfloor (u_n b + u_{n-1})/v_{n-1} \rfloor. From the properties of the floor function we have unb+un1q^vn1+vn11u_n b + u_{n-1} \leq \hat{q} v_{n-1} + v_{n-1} - 1 and therefore q^vn1unb+un1vn1+1\hat{q} v_{n-1} \geq u_n b + u_{n-1} - v_{n-1} + 1. We then get

uq^vuq^vn1bn1unbn++u0(unb+un1vn1+1)bn1=un2bn2++u0bn1+vn1bn1<vn1bn1v.\begin{aligned} u - \hat{q} v &\leq u - \hat{q} v_{n-1} b^{n-1} \\ &\leq u_n b^n + \cdots + u_0 - (u_n b + u_{n-1} - v_{n-1} + 1) b^{n-1} \\ &= u_{n-2} b^{n-2} + \cdots + u_0 - b^{n-1} + v_{n-1} b^{n-1} < v_{n-1} b^{n-1} \leq v. \end{aligned}

So uq^v<vu - \hat{q} v < v and since 0uqv<v0 \leq u - q v < v we must have qq^q \leq \hat{q}. ◼

If uu and vv are scaled appropriately, q^\hat{q} will never be too large, either.

Theorem 2. With q^\hat{q} as defined above and vn1b/2v_{n-1} \geq \lfloor b/2 \rfloor, we have q^q+2\hat{q} \leq q+2.

Proof. Assume that q^q+3\hat{q} \geq q+3. We get

q^unbun1vn1=unbnun1bn1vn1bn1uvn1bn1<uvbn1,\hat{q} \leq \frac{u_n b u_{n-1}}{v_{n-1}} = \frac{u_n b^n u_{n-1} b^{n-1}}{v_{n-1} b^{n-1}} \leq \frac{u}{v_{n-1} b^{n-1}} < \frac{u}{v - b^{n-1}},

since v=vn1bn1++v0vn1bn1+bn1v = v_{n-1} b^{n-1} + \cdots + v_0 \leq v_{n-1} b^{n-1} + b^{n-1}. We cannot have v=bn1v = b^{n-1} since that would imply q^=q=un\hat{q} = q = u_n. The relation q=u/vq = \lfloor u/v \rfloor implies q>u/v1q > u/v - 1, from which we get

3q^q<uvbn1uv+1=uv(bn1vbn1)+1.3 \leq \hat{q} - q < \frac{u}{v - b^{n-1}} - \frac{u}{v} + 1 = \frac{u}{v} \left( \frac{b^{n-1}}{v - b^{n-1}} \right) + 1.

We then have

uv2(vbn1bn1)2(vn11),\frac{u}{v} \geq 2 \left( \frac{v - b^{n-1}}{b^{n-1}} \right) \geq 2(v_{n-1} - 1),

and finally

b4q^3q=u/v2(vn11),b-4 \geq \hat{q}-3 \geq q = \lfloor u/v \rfloor \geq 2(v_{n-1}-1),

which implies vn1<b/2v_{n-1} < \lfloor b/2 \rfloor. ◼

We would expect to come even closer if we consider the 'second order' approximate quotient,

unb2+un1b+un2vn1b+vn2,\left\lfloor \frac{u_n b^2 + u_{n-1} b + u_{n-2}}{v_{n-1} b + v_{n-2}} \right\rfloor,

but how much closer? Given some approximate quotient q^\hat{q}, let us compute the corresponding second order residual

unb2+un1b+un2q^(vn1b+vn2)=r^b+un2q^vn2,u_n b^2 + u_{n-1} b + u_{n-2} - \hat{q} (v_{n-1} b + v_{n-2}) = \hat{r} b + u_{n-2} - \hat{q} v_{n-2},

where r^\hat{r} is the first order residual,

r^=unb+un1q^vn1.\hat{r} = u_n b + u_{n-1} - \hat{q} v_{n-1}.

By studying the sign of the second order residual we can now get closer to the true quotient.

Theorem 3. Let q^\hat{q} be any approximate quotient and r^\hat{r} the corresponding first order residual. Now if q^vn2>br^+un2\hat{q} v_{n-2} > b \hat{r} + u_{n-2} then q<q^q < \hat{q}.

Proof. Assume q^vn2>br^+un2\hat{q} v_{n-2} > b \hat{r} + u_{n-2}, equivalent to r^b+un2q^vn2+10\hat{r} b + u_{n-2} - \hat{q} v_{n-2} + 1 \leq 0. We then have

uq^vuq^vn1bn1q^vn2bn2=bn1(unb+un1q^vn1)+un2bn2++u0q^vn2bn2<bn1r^+un2bn2+bn2q^vn2bn2=bn2(r^b+un2q^vn2+1)0.\begin{aligned} u - \hat{q} v &\leq u - \hat{q} v_{n-1} b^{n-1} - \hat{q} v_{n-2} b^{n-2} \\ &= b^{n-1} (u_n b + u_{n-1} - \hat{q} v_{n-1}) + u_{n-2} b^{n-2} + \cdots + u_0 - \hat{q} v_{n-2} b^{n-2} \\ &< b^{n-1} \hat{r} + u_{n-2} b^{n-2} + b^{n-2} - \hat{q} v_{n-2} b^{n-2} \\ &= b^{n-2} (\hat{r} b + u_{n-2} - \hat{q} v_{n-2} + 1) \leq 0. \end{aligned}

So uq^v<0uqvu - \hat{q} v < 0 \leq u - q v which implies q<q^q < \hat{q}. ◼

Theorem 4. Let q^\hat{q} be any approximate quotient and r^\hat{r} the corresponding first order residual. Now if q^vn2br^+un2\hat{q} v_{n-2} \leq b \hat{r} + u_{n-2} then q^q+1\hat{q} \leq q+1.

Proof. Let q^vn2br^+un2\hat{q} v_{n-2} \leq b \hat{r} + u_{n-2} and assume q^q+2\hat{q} \geq q+2. Now since uqv<vu - q v < v we get

u<(q+1)v(q^1)v<q^(vn1bn1+vn2bn2+bn2)v<q^vn1bn1+q^vn2bn2+bn1vq^vn1bn1+(br^+un2)bn2+bn1v=unbn+un1bn1+un2bn2+bn1vunbn+un1bn1+un2bn2u.\begin{aligned} u &< (q+1) v \leq (\hat{q}-1) v < \hat{q} (v_{n-1} b^{n-1} + v_{n-2} b^{n-2} + b^{n-2}) - v \\ &< \hat{q} v_{n-1} b^{n-1} + \hat{q} v_{n-2} b^{n-2} + b^{n-1} - v \\ &\leq \hat{q} v_{n-1} b^{n-1} + (b \hat{r} + u_{n-2}) b^{n-2} + b^{n-1} - v \\ &= u_n b^n + u_{n-1} b^{n-1} + u_{n-2} b^{n-2} + b^{n-1} - v \\ &\leq u_n b^n + u_{n-1} b^{n-1} + u_{n-2} b^{n-2} \leq u. \end{aligned}

This claims that u<uu < u, a contradiction, so our assumption q^q+2\hat{q} \geq q+2 must have been wrong. ◼

We now have the following procedure for computing q^\hat{q}, a very close estimate to qq:

Algorithm Q. Let u=(unu1u0)bu = (u_n \ldots u_1 u_0)_b and v=(vn1v1v0)bv = (v_{n-1} \ldots v_1 v_0)_b, vn10v_{n-1} \neq 0, with n2n \geq 2 and 0u/v<b0 \leq \lfloor u/v \rfloor < b (any digit of uu can be zero and note that the only digits accessed are unu_n, un1u_{n-1}, un2u_{n-2}, vn1v_{n-1}, and vn2v_{n-2}). The algorithm computes an integer q^\hat{q} such that q^1u/vq^\hat{q}-1 \leq \lfloor u/v \rfloor \leq \hat{q} (Theorems 1 and 4).

  • Q1. Set q^(unb+un1)/vn1\hat{q} \leftarrow \lfloor (u_n b + u_{n-1})/v_{n-1} \rfloor and r^(unb+un1) mod vn1\hat{r} \leftarrow (u_n b + u_{n-1}) \text{ mod } v_{n-1}. If q^=b\hat{q} = b (division overflow when b=bTb=b_T) set q^q^1\hat{q} \leftarrow \hat{q} - 1 and r^r^+vn1\hat{r} \leftarrow \hat{r} + v_{n-1} (dealing with division overflow can be avoided by setting q^b1\hat{q} \leftarrow b-1 and r^un+un1\hat{r} \leftarrow u_n + u_{n-1} if vn1=unv_{n-1} = u_n).
  • Q2. While r^<b\hat{r} < b and q^vn2>br^+un2\hat{q} v_{n-2} > b \hat{r} + u_{n-2}, set q^q^1\hat{q} \leftarrow \hat{q} - 1 and r^r^+vn1\hat{r} \leftarrow \hat{r} + v_{n-1} (Theorem 2 assures that this while-loop is executed at most two times if vn1b/2v_{n-1} \geq \lfloor b/2 \rfloor. The check r^<b\hat{r} < b is not necessary but makes sure that we don't deal with numbers that are b2b^2 or larger in the subsequent comparison).

We can now combine Algorithm G with the just obtained knowledge of approximating the quotient in the following algorithm for long division:

Algorithm L. Given u=(um1u1u0)bu = (u_{m-1} \ldots u_1 u_0)_b, um10u_{m-1} \neq 0 and v=(vn1v1v0)bv = (v_{n-1} \ldots v_1 v_0)_b, vn10v_{n-1} \neq 0, with mn>1m \geq n > 1, this algorithm computes the quotient q=(qmnq1q0)b=u/vq = (q_{m-n} \ldots q_1 q_0)_b = \lfloor u/v \rfloor (we may have qmn=0q_{m-n} = 0 in which case qmn10q_{m-n-1} \neq 0 if m>nm > n) and the remainder rr such that u=qv+ru = q v + r, 0r<v0 \leq r < v.

  • L1. Set vdvv \leftarrow d \cdot v such that vn1b/2v_{n-1} \geq \lfloor b/2 \rfloor (letting dd be a power of two is usually the best choice). Similarly, set (umu1u0)bdu(u_m \ldots u_1 u_0)_b \leftarrow d \cdot u (ensure uu gets n+1n+1 digits, setting um=0u_m=0 if necessary).
  • L2. Set kmnk \leftarrow m - n.
  • L3. Find q^\hat{q} such that q^1(uk+nuk+1uk)b/vq^\hat{q}-1 \leq \lfloor (u_{k+n} \ldots u_{k+1} u_k)_b /v \rfloor \leq \hat{q} (use Algorithm Q described above).
  • L4. Make the update (uk+nuk+1uk)b(uk+nuk+1uk)bq^v(u_{k+n} \ldots u_{k+1} u_k)_b \leftarrow (u_{k+n} \ldots u_{k+1} u_k)_b - \hat{q} v.
  • L5. If the subtraction of step L4 produces a borrow (the result is negative) do q^q^1\hat{q} \leftarrow \hat{q} - 1 and (uk+nuk+1uk)b(uk+nuk+1uk)b+v(u_{k+n} \ldots u_{k+1} u_k)_b \leftarrow (u_{k+n} \ldots u_{k+1} u_k)_b + v.
  • L6. Set qk=q^q_k = \hat{q}.
  • L7. If k=0k=0 set ru/dr \leftarrow u/d and exit. Otherwise set kk1k \leftarrow k-1 and go to step L3.

The normalization in step L1 such that vn1b/2v_{n-1} \geq \lfloor b/2 \rfloor does two things. Firstly, it makes sure that the while-loop of the q^\hat{q}-computation executes at most two times. Secondly, the probability that the adding back in step L5 must be executed is of order 2/b2/b (a proof can be found in Knuth's book).