11 Big Integers
11.6 Division and Remainders¶
We now come to integer division, the last of the four basic arithmetic operations, and the one that is most difficult to implement efficiently. At first, it may not be apparent why division should be more difficult than, say, multiplication. After all, if we are dividing an integer \(u\) by another integer \(v\), all we are asking is “how often does \(v\) fit into \(u\)?” and this question can be easily answered by repeatedly subtracting \(v\) and counting how many iterations it takes until the result is either 0, in which case \(v\) divides \(u\) evenly, or until the result is negative, in which case \(v\) divides \(u\) with a remainder. Unfortunately, this naive “algorithm” is painfully slow — to compute \(1\,000\,000/1\), for instance, it requires a million iterations, even though the result has only seven decimal digits. What we are looking for is an algorithm that finds each digit of the result in a constant amount of time.
Before we go into the details, let’s first clarify the definition of integer division. If \(x\) and \(y\) are two positive integers, their quotient \(q=\idiv {x}{y}\) and remainder \(r=x\bmod y\) are the unique positive integers that satisfy the following two conditions:
\(\seteqnumber{0}{11.}{9}\)\begin{equation} \label {eq:divrem} x = qy + r\qquad \text {and}\qquad r < y. \end{equation}
The first part \(x=qy+r\) tells us that integer division is the inverse of multiplication, whereas the second part \(r<y\) stipulates that the remainder must be as small as possible.
-
Example 11.5. Dividing \(x=17\) by \(y=3\) is equivalent to solving the equation \(x=qy+r\) for the quotient \(q\) and the remainder \(r\). There are multiple solutions. One is \(q=1, r=14\) since we have \(17=1\cdot 3+14\), and the others are \((2,11)\), \((3,8)\), \((4,5)\), and \((5,2)\). Only the last solution \(q=5\) and \(r=2\) satisfies the additional constraint \(r<y\), so we have \(\lfloor 17/3\rfloor =5\) and \(17\bmod 3=2\).
Our goal in the remainder of this section is to find an algorithm for integer division that computes both the quotient and the remainder of two arbitrary big integers. We start with the traditional pen-and-paper method for dividing decimal numbers and then make a series of adjustments to handle other radixes and turn the method into a proper algorithm.
-
Example 11.6. The pen-and-paper method for long division computes the digits of the quotient from left to right, starting with the most significant digit. Let’s see this method in action and divide \(x=78\,539\) by \(y=382\). The whole process is usually laid out in a tabular form like this:2
\(\seteqnumber{0}{11.}{10}\)\begin{equation*} \begin{array}{r@{\,}c@{\,}l@{\extracolsep {1cm}}l} &&{\underline {\phantom {00}\mathbf {205}}}\\ 382 &\bigl )& 78539\\ &&\underline {764\phantom {00}} & 382\cdot \mathbf {2}\le 785\\ &&\phantom {0}213 & \text {subtract, bring down next digit}\\ &&\underline {\phantom {000}0\phantom {0}} & 382\cdot \mathbf {0}\le 213\\ &&\phantom {0}2139 & \text {subtract, bring down next digit}\\ &&\underline {\phantom {0}1910} & 382\cdot \mathbf {5}\le 2139\\ &&\phantom {00}229 & \text {remainder} \end {array} \end{equation*}
The three digits \(q_2q_1q_0\) of the result \(205\) are computed sequentially. For the first digit \(q_2\), we determine how often the divisor \(382\) fits into the leading digits \(785\) of the dividend. This can be accomplished using trial and error; here, it’s easy to see that \(2\cdot 382=764\) is barely smaller than \(785\), so the leftmost digit must be \(q_2=2\). Before we compute the next digit, we remove the contribution of \(q_2\) from the dividend by subtracting \(764\) from the leading digits and then append the next digit from \(x\), in this case 3. For the next digit, we have to compute \(213/382\), which immediately gives us \(q_1=0\). We again subtract the contribution of \(q_1\cdot 382=0\) and append the next digit of \(x\).
To find the final digit \(q_0\) we have to divide \(2139\) by \(382\), which requires a bit of trial and error. Since \(3\cdot 7=21\), a reasonable first guess is \(q_0=7\), but a test multiplication shows that this estimate is too large since \(7\cdot 382=2674\). We then try \(6\cdot 382=2292\), which is still too large, and \(5\cdot 382=1910\), which finally gives us a result smaller than \(2139\). The last digit is therefore \(q_0=5\), and performing one more subtraction gives us the remainder \(229\). It is easy to verify that this result is correct: \(205\cdot 382+229=78\,539\).
It is possible to generalize this procedure to numbers in any base \(B\). Assume we want to compute the quotient \(\lfloor x/y\rfloor \) and the remainder \(x\bmod y\) of two arbitrary base-B integers \(x\) and \(y\). If \(x\ge y\), we denote by \(m\) the number of digits in the divisor \(y\) and by \(n\) the number of additional digits in the dividend \(x\); in other words, we have
\[ x=(x_{n+m-1}\ldots x_0)_B,\quad \text {and}\quad y=(y_{m-1}\ldots y_0)_B. \]
With this convention, the quotient has up to \(n+1\) digits and the remainder up to \(m\) digits.
As we will explain in more detail in Box 11.4, the first digit of the quotient \(q_{n}\) can be found by taking the \(m\) most significant digits of \(x\) and dividing them by \(y\):
\[ q_{n} = \bigl \lfloor (x_{n+m-1}\dots x_{n})_B / y\bigr \rfloor . \]
As in the pen-and-paper method, we then subtract \(q_{n}y\) from \((x_{n+m-1}\dots x_{n})_B\) to compute the remainder of the first step
\[ r_{n} = (x_{n+m-1}\dots x_{n})_B - q_{n}y. \]
We can then find the next digit of the quotient by extending \(r_{n}\) with the \(n-1\)th digit from \(x\) and again dividing the result by \(y\). This process is repeated until all \(n+1\) digits \(q_n,\ldots ,q_0\) have been determined.
The resulting algorithm Divide is shown in Algorithm 11.6. Each iteration starts by computing \(t=rB+x_i\), which appends the next digit \(x_i\) of the dividend to the remainder \(r\) of the previous step. Before the loop, \(r\) is initialized to the \(m-1\) most significant digits of \(x\).
Algorithm 11.6
-
Example 11.7. To see that Divide is just a generalized version of the normal pen-and-paper method for long division, let’s use it to divide \(x=78\,539\) by \(y=382\). The divisor \(y\) has three digits, so \(m=3\), and the dividend \(x\) has two additional digits, so \(n=2\). The algorithm now computes the digits \(q_2,q_1,q_0\) in this order. The following table shows the values of all variables used inside the for loop of Divide:
\(\seteqnumber{0}{11.}{10}\)\begin{equation*} \begin{array}{rrrrr} i & t & q_i & r \\ & & & 78 \\ 2 & 785 & 2 & 21 \\ 1 & 213 & 0 & 213\\ 0 & 2139 & 5 & 229\\ \end {array} \end{equation*}
In the first row, the variable \(r\) is initialized to the \(m-1\) most significant digits of \(x\). Each of the following steps then extends the previous value of \(r\) by appending the next digit from \(x\), computes the \(i\)th digit \(q_i=\idiv {t}{y}\), and updates the remainder to \(t-q_iy\).
Notice that algorithm Divide computes each digit of the quotient exactly, without having to resort to the usual trial and error we associate with long division. This assumes, however, that we can actually evaluate the expression \(\lfloor t/y\rfloor \) directly, which isn’t always the case — after all, \(t\) and \(y\) are still arbitrarily large big integers. Nevertheless, computing \(\idiv {t}{y}\) is significantly easier than the original problem \(\idiv {x}{y}\) since we know that the result must be a single digit; we will see why this is important shortly.
Let’s start with a special case in which \(\lfloor t/y\rfloor \) can in fact be evaluated directly. If the divisor \(y\) consists of a single digit, we have \(m=1\) and the division procedure simplifies to the Divide-Digit algorithm shown in Algorithm 11.7. The key property of Divide-Digit is that all quantities used inside the algorithm have only one or two digits: \(y\) and \(r\) are obviously single-digit numbers, and it’s not hard to show that the intermediate value \(t\) has two digits at most (see Exercise 11.10).
Algorithm 11.7
Because Divide-Digit operates on quantities with at most two digits, its implementation for big integers is straightforward. The divideByDigit() method in Listing 11.15 divides a BigNat by a single digit y and returns the quotient and the remainder. The divisor y is of type int and the remaining quantities used by the algorithm are stored in variables of type long. We can therefore use ordinary 64-bit arithmetic compute \(t\), \(q_i\) and \(r\).
Listing 11.15ch11β―/β―BigNat
// Divide 'this' by 'y', return the quotient and remainder.
private BigNat[] divideByDigit(int y) {
int[] quot = new int[size];
long r = 0;
for (int i = size - 1; i >= 0; i--) {
long t = (r << DIGIT_BITS) + (long) digit(i);
long qi = t / y;
quot[i] = (int) qi;
r = t - y * qi;
}
return new BigNat[]{new BigNat(quot), fromUnsignedLong(r)};
}
Box 11.4: Understanding Integer DivisionWhy and how does the division algorithm in Algorithm 11.6 work? It’s easy to understand its mechanics by starting with the inverse operation, the multiplication of two integers:
In the division problem we now take the product \(x=436749\) and the second factor \(y=739\) and try to recover the first factor \(q=\lfloor x/y\rfloor =591\).
Since \(x\) and \(y\) have 6 and 3 digits, respectively, their quotient \(q\) has at most three digits \(q_2q_1q_0\), and the multiplication procedure shown above makes it clear that the most significant digit \(q_2\) only affects the four highest digits of \(x\). We therefore have
\(\seteqnumber{0}{11.}{10}\)\begin{equation} \label {eq:divfirstdigit} 4367 = q_2\cdot 739 + r_2, \end{equation}
where \(r_2\) is the contribution of multiplying 739 by the remaining digits \(q_1\) and \(q_0\). It is possible to show that \(r_2<y\)3, so we can solve Eq. (11.11) for \(q_2\) and \(r_2\) using integer division:
\(\seteqnumber{0}{11.}{11}\)\begin{equation*} q_2 = \lfloor 4367 / 739\rfloor = 5,\qquad r_2 = 4367\bmod 739 = 672. \end{equation*}
We can remove now remove \(q_2\)’s contribution to the product by subtracting \(q_2\cdot 739\) from the leading digits of \(x\), which gives us \(67249\). The problem of finding the remaining digits of \(q\) is now equivalent to inverting the following multiplication:
This problem can be solved in the same way, which gives us the remaining digits \(q_1\) and \(q_0\).
In this example, the final remainder is \(r_0=0\) because 739 divides 436749 evenly. In other division problems, \(r_0\) is equal to the remainder of the original division \(r_0=x \bmod y\) — it’s the remaining term after removing the contributions of each of \(q\)’s digits to \(x\).
3 The proof similar to the one in Exercise 11.5.
Finding Digits using Trial and Error
Let’s now turn to the problem of implementing Divide when the divisor \(y\) is not a single digit but an arbitrary big integer. The main problem in this case is computing the \(i\)th digit \(q_i=\lfloor t/y\rfloor \) of the quotient. We know that the correct result \(q_i\) satisfies two conditions:
\[ q_i < B,\qquad \text {and}\qquad q_i y + r = t,\quad \text {with $r < y$}. \]
We can therefore find \(q_i\) by trying all possible digits to find the largest one that satisfies \(q_i y \le t\):
Unfortunately, this approach is prohibitively slow, especially if the radix is large: In our case, it requires up to \(B=2^{31}\) iterations to find just a single digit of the quotient.
We can drastically improve this algorithm by initializing \(q_i\) to a value that is closer to the true result \(\idiv {t}{y}\). The most common strategy for estimating \(q_i\) is to take the two leading digits of \(t\) and divide them by the leading digit of \(y\):
\(\seteqnumber{0}{11.}{11}\)\begin{equation} \label {eq:qhat} \hat q=\left \lfloor \frac {t_{m}\cdot B + t_{m-1}}{y_{m-1}}\right \rfloor . \end{equation}
It can be shown that this estimate always satisfies \(\hat q\ge \lfloor t/y\rfloor \), in other words, it is never smaller than the actual result.4 We can therefore use the same while loop as before to adjust \(q_i\) until it satisfies \(q_i y \le t\):
4 For a proof, see [60, Section 4.3.1].
If \(y_{m-1}\) is small, Eq. (11.12) can overestimate the true result by a lot, so we initialize \(q_i\) to the minimum of \(\hat q\) and \(B-1\).
If we use this technique to compute \(\lfloor t/y\rfloor \) in algorithm Divide, we obtain the updated algorithm for dividing big integers shown in Algorithm 11.8.
Algorithm 11.8
Compute the quotient \(\idiv {x}{y}\) and the remainder \(x\bmod y\) of two integers \(x=(x_{n+m-1}\ldots x_0)_B\) and \(y=(y_{m-1}\ldots y_0)_B\). This is a variation of Algorithm 11.6 that uses only two-digit divisions internally.
Of course, this algorithm is only better than simply trying all possible digits if \(\hat q\) is close to the true digit \(q_i\). This isn’t always the case, as the following example demonstrates.
-
Example 11.8. Let’s use Divide-Long to compute \(\idiv {1000}{19}\) in the decimal system. We have \(m=2\) and \(n=2\) and the algorithm performs three iterations to compute the digits of the quotient \((q_2q_1q_0)_{10}\). The following table lists the intermediate values computed:
\(\seteqnumber{0}{11.}{12}\)\begin{equation*} \begin{array}{crccc} i & t & \hat q & q_i & r\\ 2 & 010 & \lfloor 01 / 1\rfloor = 1 & 0 & 10\\ 1 & 100 & \min (\lfloor 10 / 1\rfloor , 9) = 9 & 5 & 5\\ 0 & 050 & \lfloor 05 / 1\rfloor = 5 & 2 & 12\\ \end {array} \end{equation*}
The columns for \(\hat q\) and \(q_i\) show the original estimate for \(\idiv {t}{y}\) and the true result; the difference between these two numbers is the number of correction steps performed by the inner while loop. The result is therefore \(\idiv {1000}{19}=052\) with a remainder of \(12\), and the total number of correction steps performed in this case is \(1+4+3=8\).
The algorithm performs much better for other divisors. For example, here are the intermediate values when computing \(\idiv {1000}{52}\):
\(\seteqnumber{0}{11.}{12}\)\begin{equation*} \begin{array}{crccc} i & t & \hat q & q_i & r\\ 2 & 010 & \lfloor 01 / 5\rfloor = 0 & 0 & 10\\ 1 & 100 & \lfloor 10 / 5\rfloor = 2 & 1 & 48\\ 0 & 480 & \lfloor 48 / 5\rfloor = 9 & 9 & 12\\ \end {array} \end{equation*}
In this case, the estimate \(\hat q\) is always close to the true digit \(q_i\) and only a single correction step is necessary in total.
The reason why Divide-Long has to perform more corrections when dividing by 19 than when dividing by 51 is that we use only the leftmost digit \(y_{m-1}\) of the divisor to compute \(\hat q\). If \(y_{m-1}\) is large compared to the remaining digits (as for \(y=52\)), omitting the remaining digits of \(y\) introduces only a small error and \(\hat q\) is close to the true digit. But if \(y_{m-1}\) is relatively small (as for \(y=19\)), we introduce a large error and the estimate \(\hat q\) will be off. In the decimal system, dividing 19, 199, etc. is the worst case scenario that requires up four correction steps per digit. Similarly, the worst case in base-\(B\) is dividing by numbers of the form \((1 (B-1) \ldots (B-1))_B\), in which case up to \(\lfloor B/2\rfloor \) corrections are necessary — that’s still one billion iterations per digit when dividing certain big integers!
Fortunately, there is a simple trick that can be used to ensure that we never need to apply more than a few correction steps. The idea is to multiply both the dividend and divisor of the original division \(x/y\) by a common factor \(\alpha \). This doesn’t change the result of the division since
\(\seteqnumber{0}{11.}{12}\)\begin{equation*} \frac {\alpha x}{\alpha y} = \frac {x}{y}, \end{equation*}
but if we choose \(\alpha \) in such a way that the leading digit of \(\alpha y\) is sufficiently large, we can avoid the worst-case behavior of algorithm Divide-Long. It can be shown that the choice
\(\seteqnumber{0}{11.}{12}\)\begin{equation} \label {eq:divscaling} \alpha = \left \lfloor \frac {B-1}{y_{m-1}+1}\right \rfloor \end{equation}
guarantees that the leading digit of \(\alpha y\) is greater or equal than \(\idiv {B}{2}\) and that in this case no more than two correction steps are required, regardless of the radix \(B\). This is a remarkable improvement that only requires two additional multiplications to compute \(\alpha x\) and \(\alpha y\) and one single-digit division to recover the remainder of \(x\) and \(y\) from the remainder of \(\alpha x\) and \(\alpha y\).
Implementing Long Division
We can finally tie everything together and finish the implementation of integer division. We start with a top-level method for dividing two arbitrary big integers \(x\) and \(y\). We distinguish the following four cases:
-
1. If \(y=0\) we raise an error since division by zero is undefined.
-
2. If \(y\) has a larger magnitude than \(x\), the quotient \(\idiv {x}{y}\) is always 0 and the remainder is \(x\bmod y=x\).
-
3. If \(y\) has just a single digit we use algorithm Divide-Digit.
-
4. If \(y\) has two or more digits we use algorithm Divide-Long.
The divideAndRemainder() function in Listing 11.16 handles division by zero and division by a larger divisor directly and delegates the other two cases to divideByDigit() and divideLong(). All division functions return an array of two BigNats containing the quotient \(\idiv {x}{y}\) and the remainder \(x\bmod y\). For convenience, we also define two methods dividedBy() and remainder() that call divideAndRemainder() and return just the quotient or just the remainder.
Listing 11.16ch11β―/β―BigNat
// Compute the quotient and remainder of this and 'y'.
public BigNat[] divideAndRemainder(BigNat y) {
if (y.equals(ZERO)) {
throw new ArithmeticException("Division by zero");
} else if (this.compareTo(y) < 0) {
return new BigNat[]{ZERO, this};
} else if (y.size == 1) {
return divideByDigit(y.digit(0));
} else {
return divideLong(y);
}
}
public BigNat dividedBy(BigNat y) {
return divideAndRemainder(y)[0];
}
public BigNat remainder(BigNat y) {
return divideAndRemainder(y)[1];
}
The implementation of the long division algorithm in Listing 11.17 consists of three methods. The top-level method is divideLong(), which scales the dividend and divisor according to Eq. (11.13) and then calls divideLongInternal() to perform the actual division. Since the scaling factors cancel out, the resulting quotient can be returned as is, but the remainder must be scaled back by dividing it by the scaling factor.
Listing 11.17ch11β―/β―BigNat
private BigNat[] divideLong(BigNat y) {
int f = (int) ((BASE - 1) / (y.digit(y.size - 1) + 1));
BigNat[] divRem = this.timesDigit(f)
.divideLongInternal(y.timesDigit(f));
return new BigNat[]{divRem[0],
divRem[1].divideByDigit(f)[0]};
}
private BigNat[] divideLongInternal(BigNat y) {
int n = size - y.size;
int m = y.size;
var r = slice(n + 1, m - 1);
int[] quot = new int[n + 1];
for (int i = n; i >= 0; i--) {
var t = r.appendDigit(digit(i));
long head = ((long) t.digitOr0(m) << DIGIT_BITS)
+ t.digitOr0(m - 1);
long qi = Math.min(head / y.digit(m - 1), BASE - 1);
BigNat yTimesQ = y.timesDigit((int) qi);
while (t.compareTo(yTimesQ) < 0) {
qi--;
yTimesQ = yTimesQ.minus(y);
}
quot[i] = (int) qi;
r = t.minus(yTimesQ);
}
return new BigNat[]{new BigNat(quot), r};
}
private BigNat appendDigit(int digit) {
int[] shifted = new int[size + 1];
System.arraycopy(digits, offset, shifted, 1, size);
shifted[0] = digit;
return new BigNat(shifted);
}
JDK: Math, System
•BigNat: BigNat(), compareTo(), digit(), digitOr0(), divideByDigit(), minus(), slice(), timesDigit()
The divideLongInternal() method implements the Divide-Long algorithm in Algorithm 11.8. We first initialize the variables n, m, and r to the same values they have in Divide-Long and then allocate an array quot that holds the digits of the quotient. Each iteration of the for loop then computes one of these digits.
Inside the loop, we use two simple optimizations to speed up our implementation. First, we note that the expression \(t = rB + x_i\) in the original algorithm shifts the digits of \(r\) one place to the left and inserts the next digit \(x_i\) of the dividend as the least significant digit. Instead of performing a multiplication followed by an addition, the appendDigit() method implements this operation by copying the digits to a new array. A second optimization is used in the while loop that adjusts \(q_i\), where we compute the product \(q_i y\) just once before the loop and then adjust it by subtracting \(y\) every time we decrement \(q_i\).
The implementation of long division can be made more efficient by ensuring that no new big integers or arrays of digits are allocated inside the main loop. We can do this by creating two arrays before the loop that hold the digits of t and yTimesQ and then rewriting methods such appendDigit(), compareTo(), and minus() to operate on those arrays.
Exercises
Exercise 11.9.The Divide algorithm on page (Page for Algorithm 11.6) requires us to compute \(\lfloor t/y\rfloor \) to find each digit of the quotient. Is it possible to perform this computation by calling Divide recursively?
Exercise 11.10. Show that in each iteration of Divide-Digit (Algorithm 11.7) the variable \(t\) satisfies \(t<B^2\) and therefore has at most two digits. Hint: Use the fact that we have \(r=t\bmod y\) and therefore \(r<y\le B-1\).
Exercise 11.11. Manually evaluate Divide-Long\((785, 18)\) in the decimal system. Repeat the computation with both arguments scaled by the factor in Eq. (11.13) and compare the total number of correction steps that are necessary.