11 Big Integers
11.9 Problems¶
Analysis of Karatsuba Multiplication
\(\star \) Exercise 11.17. The runtime complexity of many divide-and-conquer algorithms can be analyzed using the so-called master theorem. Assume we are given a problem of size \(n\) (in the case of multiplication, \(n\) is the number of digits of the two numbers being multiplied) and let \(T(n)\) denote the time required to solve this problem using a particular algorithm. If the algorithm recursively splits the problem into \(a\) smaller problems of size \(n/b\), we can write \(T(n)\) as
\(\seteqnumber{0}{11.}{17}\)\begin{equation} \label {eq:masterbase} T(n) = aT\left (\frac {n}{b}\right ) + \left (\begin{array}{@{}c@{}} \text {Time required for splitting problem}\\ \text {and combining solutions} \end {array}\right ). \end{equation}
The master theorem now states that \(T(n)\) has the asymptotic complexity
\(\seteqnumber{0}{11.}{18}\)\begin{equation} T(n) = O(n^{\log _b a}), \end{equation}
provided the second term on the right-hand side of Eq. (11.18) is sufficiently small.5
Here is a simple example. In Section 3.1 we discussed a simple recursive algorithm for finding the smallest element in an array of integers. The algorithm splits the array into two parts, recursively finds the minimum of each, and then returns the smaller of the two values. If we split the array into two equal parts, we have \(a=2\) and \(b=2\), so
\[ T(n) = 2\cdot T(n / 2) + C, \]
where \(C\) is a constant that measures the work done before and after the recursive calls, which is independent of \(n\). According to the master theorem, the complexity of the algorithm is therefore
\[ T(n) = O(n^{\log _2 2}) = O(n). \]
-
a) What are the values of \(a\) and \(b\) for the recursive multiplication scheme in Eqs. (11.7) and (11.8)? Assume that master theorem can be applied. What is the resulting asymptotic complexity of \(T(n)\)?
-
b) What are the values of \(a\) and \(b\) for Karatsuba’s algorithm? Again, assume that master theorem can be applied and derive the asymptotic complexity of \(T(n)\).
5 For the technical details and a more general version of the master theorem, see [25, Chapter 4].
Integer Powers
Exercise 11.18.Consider the problem of computing \(x^m\) for two integers \(x\) and \(m\). The obvious algorithm computes the product \(x\cdot x\cdot \dots \cdot x\), where \(x\) is repeated \(m\) times, which requires a total number of \(m-1\) multiplications. A much faster algorithm for computing integer powers can be obtained by repeated squaring. It’s easiest to see how it works by considering an example. If we want to compute \(x^{19}\), we start with the binary expansion of the exponent \(19=10011_2=1+2+16\) and then write \(x^{19}\) as a product of three terms of the form \(x^e\), where \(e\) is a power of 2:
\(\seteqnumber{0}{11.}{19}\)\begin{equation*} x^{19}=x^{1+2+16}=xx^2x^{16}. \end{equation*}
Each of the terms \(x^e\) can be obtained by repeatedly squaring the previous term: \(x^2\) is obtained by squaring \(x\) once, and \(x^{16}\) is obtained by squaring \(x^2\) three times:
\(\seteqnumber{0}{11.}{19}\)\begin{equation*} x\xrightarrow {\text {square once}} x^2 \xrightarrow {\text {square 3 times}} x^{16}. \end{equation*}
This method reduces the problem of computing \(x^{19}\) from 18 multiplications to only six multiplications: four multiplications to compute the terms \(x^e\) and two multiplications to accumulate their product.
We can turn this idea into a proper algorithm as follows. Assume we want to compute \(x^m\) for some integer \(x\) and and a positive exponent \(m\). As in the example we just discussed, the first step is to write \(x^m\) as a product of smaller powers of \(x\) using the binary representation of the exponent
\(\seteqnumber{0}{11.}{19}\)\begin{equation*} m = (m_{M-1} \dots m_1m_0)_2 = m_{M-1} 2^{M-1} + \dots + m_1 2^1 + m_0 2^0. \end{equation*}
Using this expansion, we can use the fact that \(x^{a+b}=x^ax^b\) to rewrite \(x^m\) as a product of smaller powers:
\(\seteqnumber{0}{11.}{19}\)\begin{align*} x^m &= x^{m_{M-1} 2^{M-1} + \dots + m_1 2^1 + m_0 2^0}\\ &=x^{m_{M-1} 2^{M-1}} \cdots x^{m_1 2^1} x^{m_0 2^0}; \end{align*} This expression can be simplified further since each bit \(m_k\) is either 0 or 1, which means that each term can be written as
\(\seteqnumber{0}{11.}{19}\)\begin{equation*} x^{m_k 2^k} = \begin{cases} x^{2^k} & \text {if $m_k=1$};\\ 1 & \text {otherwise}. \end {cases} \end{equation*}
Only the terms for which \(m_k=1\) contribute to the final product, so we have our final result:
\(\seteqnumber{0}{11.}{19}\)\begin{equation} x^m = \prod _{0\le k<M, m_k=1} x^{2^k}. \end{equation}
The \(\prod {}\) operator is essentially the mathematical notation for a loop that computes a product of the terms on its right. To translate this expression into an actual loop, we iterate over the nonzero binary digits of the exponent \(m\), and for each such digit \(k\) multiply the result by the term \(x^{2^k}\). The successive terms \(x^{2^k}\) can be easily computed as we go, starting with \(x\) as the first term and computing the following terms \(x^2, x^4, \dots \) by repeated squaring.
-
b) Use this idea to implement an algorithm that computes \(x^k\) for two integers \(x\) and \(k\).
-
c) How many multiplications does this algorithm perform when computing \(x^k\)?
This algorithm for computing \(x^k\) makes only two assumptions: the exponent \(k\) is a nonnegative integer and the base \(x\) is a mathematical object that can be multiplied. We can therefore use the same approach to compute integer powers of other mathematical objects, including real and complex numbers, fractions, and matrices.
-
d) Implement a method BigNat.pow() that computes the integer power of a big integer. Assume that the exponent is of type int.
Decimal Digits
Exercise 11.19. A positive integer with decimal digits \(abcde\dots \) is called polydivisible if every number formed by the first \(n\) digits is divisible by \(n\). For example, the number 62 485 is polydivisible because the first digit 6 is divisible by 1, the first two digits 62 are divisible by 2, and the numbers 624, 6248, and 62 485 are divisible by 3, 4, and 5, respectively.
-
a) There are nine polydivisible numbers with a single digit: 1, 2, 3, 4, 5, 6, 7, 8, and 9. How many polydivisible numbers with two digits are there?
-
b) Write a program that computes all polydivisible numbers. How many polydivisible numbers are there? How many digits does the largest polydivisible number have? What is the first polydivisible number that contains all decimal digits exactly once?
Exercise 11.20.A rational number is a quotient of the form \(n/m\), where \(n\) and \(m\) are integers. Every rational number has a decimal representation that consists of an integer part and a fractional part. For instance, we have \(7/4=1.75\) with the integer part 1 and the fractional part \(.75\).
-
a) Given two positive integers \(n\) and \(m\), devise an algorithm that computes a certain number of decimal digits of \(n/m\). For example, computing 5 decimal digits of \(7/3\) should produce 2.3333.
-
b) Use this algorithm to compute the first few thousand decimal digits of \(1/998\,001\).
Rational numbers have the following property: their decimal expansion is either finite or eventually produces a sequence of digits that repeats indefinitely; the length of this repeating sequence of digits is called the number’s period. For \(1/4\) the decimal expansion is \(0.25\), which is finite. For \(25/12\), however, the decimal expansion \(2.08\overline {3}\) continues forever; since the repeating sequence consists of a single digit, its period is 1.
-
c) Extend your algorithm to find the shortest sequence of repeating digits in the decimal expansion of \(n/m\). What is the period of \(987\,654\,321/123\,456\,789\)? Hint: Don’t search for repetitions in the sequence of decimal digits. Instead, determine when the algorithm returns to a previous state.
Durations
Exercise 11.21.Time durations are traditionally specified in a way that resembles a positional number system. A duration like “three weeks, two days, one hour, five minutes, and ten seconds” can be represented as a sequence of “digits” \(3, 2, 1, 5, 10\) in which each position represents one unit of time. Compared to ordinary number systems, the main difference is that each digit now has its own radix: Seconds and minutes are stored in base-60, hours in base-24, days in base-7, and assuming that weeks are the largest unit we use, we can choose as its radix any sufficiently large integer. Number systems in which each position has its own radix are also known as mixed-radix systems. In this exercise, we will see how arithmetic works in such systems.
-
a) Design a class that stores durations consisting of weeks, days, hours, minutes, and seconds. It should be possible to represent both positive and negative durations. The format should be normalized, so that two durations only have the same representation if they also represent the same time span.
Note: Since the focus of this exercise is on algorithms for mixed-radix systems, durations should be stored as an array of “digits”, not as a single number such as the total number of seconds. For instance, a duration of 70 seconds should be stored as 1 minute and 10 seconds.
-
b) Explain how to compare two durations for equality and how to implement the Comparable interface.
-
c) Define factory methods that construct durations representing a certain number of weeks, days, hours, minutes, or seconds. An example of such a factory method would be a function days() that is declared as follows:
public static Duration days(int n);
In this case, calling days(-8) should return a duration that represents minus 1 week and 1 day.
Also define corresponding methods for querying the number of weeks, days, etc. stored in a duration.
To add and subtract durations, we can use algorithms similar to those for big integers. The main difference is that each unit of time has its own radix, so subtracting 2 seconds from 1 day produces 23 hours, 59 minutes, and 58 seconds.
-
d) Implement methods for adding and subtracting durations.
In our discussion so far we have represented durations using a mixture of units such as “1 hour and 30 minutes”, but it is also possible to express the same durations using single units of time such as “1.5 hours” or “90 minutes” or “5400 seconds”.
-
e) Implement methods such as totalWeeks() or totalDays() that convert a duration to the given unit of time, returning the result as a double.
Computing Pi
Exercise 11.22. The problem of computing numerical approximations to \(\pi \) — the ratio of a circle’s circumference to its diameter — has a long and storied history. The first real algorithm for estimating \(\pi \) is due to the Greek mathematician Archimedes of Syracuse.
The idea behind Archimedes’ algorithm is to approximate circles using many-sided polygons. For example, given a circle of radius 1, we can construct two regular \(n\)-sided polygons, one inscribed in the circle and one circumscribed around it, as illustrated in Fig. 11.1. It is clear that the circle’s circumference, which is \(2\pi \), must larger than the perimeter of the inscribed polygon, which we call \(p_n\), but smaller than the perimeter of the circumscribed polygon, which we call \(P_n\):
\(\seteqnumber{0}{11.}{20}\)\begin{equation*} p_n < 2\pi < P_n. \end{equation*}
The more sides the polygons have, the closer they fit the circle, so as we increase \(n\), the values \(p_n\) and \(P_n\) converge toward the true circumference \(2\pi \).
Archimedes’ crucial insight was that it is possible to compute \(p_n\) and \(P_n\) for large values of \(n\) without actually constructing the corresponding polygons. Instead, we can use the following formulas to repeatedly double the number of sides:
\(\seteqnumber{0}{11.}{20}\)\begin{equation} \label {eq:archimedes} p_{2n} = \frac {2p_nP_n}{p_n + P_n},\qquad P_{2n} = \sqrt {p_{2n} P_n}. \end{equation}
This formula is also known as Archimedes’ recurrence formula (see Box 11.5 for a proof). The first identity states that \(p_{2n}\) is the harmonic mean of \(p_n\) and \(P_n\), and the second identity states that \(P_{2n}\) is the geometric mean of \(p_{2n}\) and \(P_n\).
For small \(n\) we can derive the values of \(p_n\) and \(P_n\) geometrically. For \(n=4\) the inscribed square has a side length of \(1/\sqrt {2}\) and the circumscribed square has a side length of 1, so we have \(p_4 = 2\sqrt {2}\) and \(P_4= 4\). Similarly, for \(n=6\) the inscribed and circumscribed hexagons have a side length of \(1/2\) and \(1/\sqrt {3}\), respectively, so the perimeters are \(p_6 = 3\) and \(P_6 =2\sqrt {3}\). Using these starting values we then repeatedly apply Archimedes’ recurrence formula to compute increasingly good estimates of \(\pi \).
-
a) Write a program that computes an approximation to \(\pi \) using Archimedes recurrence formula. Estimate the number of correct decimal digits this algorithm produces per iteration.
Box 11.5: Derivation of Archimedes’ AlgorithmArchimedes’ recurrence formula in Eq. (11.21) is easy to prove using basic trigonometry. The idea is to take a circle of radius 1 and construct inscribed and circumscribed \(n\)-gons as indicated in the following image:
Each segment consists of two right triangles, so the perimeter \(p_n\) of the inscribed polygon and the perimeter \(P_n\) of the circumscribed polygon can be expressed as
\(\seteqnumber{0}{11.}{21}\)\begin{equation*} P_n = 2n\tan \pi /n,\qquad p_n = 2n\sin \pi /n. \end{equation*}
The formulas for \(P_{2n}\) and \(p_{2n}\) now follow from the following trigonometric identities for half angles:
\(\seteqnumber{0}{11.}{21}\)\begin{align*} \tan (\alpha /2)=\frac {\sin \alpha }{1+\cos \alpha } = \frac {1-\cos \alpha }{\sin \alpha },\qquad \sin (\alpha /2)= \sqrt {\frac {1-\cos \alpha }{2}}. \end{align*} By substituting \(\alpha =\pi /n\) we get:
\(\seteqnumber{0}{11.}{21}\)\begin{align*} P_{2n} &= 4n\tan (\alpha /2) = 4n \frac {\sin \alpha }{1 + \cos \alpha } = 4n \frac {\sin \alpha \tan \alpha }{\sin \alpha +\tan \alpha } = \frac {2P_np_n}{P_n + p_n}, \intertext {and} p_{2n} &= 4n\sin (\alpha /2) = \sqrt {8n^2(1-\cos \alpha )} = \sqrt {8n^2\sin \alpha \tan (\alpha /2)} = \sqrt {P_{2n}p_n}. \end{align*}
Exercise 11.23. For many centuries, Archimedes’ method remained the best known method for computing \(\pi \), but its slow convergence and the amount of computational work it requires made it hard to use in practice. It wasn’t until the development of calculus in the 17th century — almost two millennia later — that more efficient methods for approximating \(\pi \) finally became available.
The most important formula from this era is Machin’s formula, which relates \(\pi \) to the arctan function:
\(\seteqnumber{0}{11.}{21}\)\begin{equation} \label {eq:machin} \frac {\pi }{4}=4\arctan \frac {1}{5}-\arctan \frac {1}{239}. \end{equation}
The arctan function itself can be approximated very efficiently using its Taylor expansion:
\(\seteqnumber{0}{11.}{22}\)\begin{equation} \label {eq:arctan} \arctan \left (\frac {1}{x}\right ) = \frac {1}{x}-\frac {1}{3x^3} + \frac {1}{5x^5}-\cdots . \end{equation}
For \(x>1\) this series expansion converges very quickly since the denominators grow rapidly.
-
a) Use Machin’s formula to manually compute an approximation to \(\pi \); use the first two terms of the expansion in Eq. (11.23).
It is easy to adapt Machin’s formula to compute digits of \(\pi \) using integer arithmetic only. First, we scale Eq. (11.22) by multiplying both sides by a large integer \(U\):
\(\seteqnumber{0}{11.}{23}\)\begin{equation} \label {eq:machin-scaled} \pi U=16U\arctan \frac {1}{5}-4U\arctan \frac {1}{239}. \end{equation}
If we set \(U=10^k\), the integer part of the right-hand side contains the first \(k\) decimal digits of \(\pi \). To compute this integer part, we first note that \(A\cdot \arctan 1/x\) is
\(\seteqnumber{0}{11.}{24}\)\begin{equation*} A\arctan \left (\frac {1}{x}\right ) = \frac {A}{x}-\frac {A}{3x^3} + \frac {A}{5x^5}-\cdots . \end{equation*}
Now, ignoring the fractional parts of all terms, we replace the divisions on the right-hand side with integer divisions:
\(\seteqnumber{0}{11.}{24}\)\begin{equation} \label {eq:arctan1x} A\cdot \arctan \left (\frac {1}{x}\right ) \approx \left \lfloor \frac {A}{x}\right \rfloor - \left \lfloor \frac {A}{3x^3}\right \rfloor + \left \lfloor \frac {A}{5x^5}\right \rfloor - \cdots \end{equation}
Since the denominators \((2k + 1)x^{2k + 1}\) grow very quickly, almost all of the resulting terms are zero, so we can truncate the series after a finite number of terms.
-
b) Use this approximation to manually compute \(100\pi \).
-
c) Implement Eq. (11.24) and Eq. (11.25) to obtain an algorithm that estimates \(U\pi \) for a given value of \(U\).
In general, we have to be careful when ignoring the fractional parts in a computation like this, but in the case of Machin’s formula, the rounding errors turn out to be fairly minor.
-
d) Use your implementation of Machin’s algorithm to compute \(10^{10000}\pi \). Find a table of \(\pi \) online; how many digits of your result are correct?