charlesreid1.com research & teaching

Project Euler Problem 1

Posted in Mathematics

permalink

Table of Contents

  • Overview: The Problem
  • Why This Problem?
  • Going Deeper: An Example

Overview: The Problem

Project Euler is a website that provides mathematically-oriented programming problems. There are many (over 500) and they are a rich source of profound mathematical insights.

I have been considering a writeup that goes deep into a particular problem, so why not do it with problem 1?

Problem 1 of Project Euler asks:

Find the sum of all the multiples of 3 or 5 below 1000.

It is a pretty simple task - one of the first things covered in a decent programming course is the assortment of mathematical operators, including the modular operator and multiplication operator, useful here.

This problem is also a familiar problem in another guise - any computer science student who has solved a variant of the fizz buzz problem will recognize this task (in the fizz buzz problem, you print fizz every time a number is divisible by 3 and buzz every time it is divisible by 5, etc.)

It is a deceptively simple problem. In fact, it is so easy to solve with a computer that you almost lose a sense of what the manual process would look like. How might we perform this task by hand?

Finally, it is an example of a problem in which we are trying to find the number of outcomes of several classes of events, and some of the events are labeled with both classes. This means it will be important to learn and apply the Inclusion-Exclusion Principle. (Fortunately, this principle is fairly straightforward to apply.)

We will get to the algorithm to counting these factors by hand, and handling more complicated constraints as well, but first I'll address why this problem - this task - was deemed important enough to be the very first step that nearly everyone takes on their epic (or... not so epic) Project Euler journey.

Why This Problem?

The central task in this problem is to find multiples of a number \(k\), and count them. The task is simple enough (unlike later Project Euler questions, which can be downright frightening at times), but still requires knowledge of loops, operators, and basic algorithm design. It's no litmus test for whether you can solve Problem 100, but it gets you started.

The task at the heart of this problem - iterating through a list of multiples of a number - is at the heart of the Sieve of Eratosthenes algorithm, which in turn is at the heart of applied number theory. While it may not be the algorithm applied in practice, it is the first and most important algorithm number theorists learn.

It's also a first lesson in the subtleties of Project Euler problems - the eager but naive algorithm designer will count all multiples of 3, then all multiples of 5, forgetting that some repeat.

Project Euler Fail

Welcome to Project Euler.

Going Deeper: An Example

It's true that this problem seems a bit boring on its face. But let's dive deeper. Suppose I asked you to find the number of multiples of the integers 3 and 4, but not the integer 5, below 2,001 - and to do so without explicitly enumerating them with a computer.

To do this, we can express the problem in set notation. We have three sets, A, B, C, containing multiples of 3, 4, and 5, respectively. In set theory language, we wish to find

$$ ( A \bigcup B ) \backslash C $$

We can start by counting the sets A and B, as well as accounting for \(A \bigcap B\) (numbers that are multiples of both a and b).

Next, we can count \(A \bigcap C\) and \(B \bigcap C\), which are the multiples of a and b that we counted that we should not have because they are multiples of c.

Finally, we cannot forget \(A \bigcap B \bigcap C\) - numbers that have a, b, and c as multiples. This case is a bit tricky. Any item that is in \(A \bigcap B \bigcap C\) has already been removed - twice. The first time was when it was removed because it was in \(A \bigcap C\), and the second time was when it was removed because it was in \(B \bigcap C\). Therefore, we must add each of these items back in, to account for the double-removal and ensure these items are only removed once.

So we will add the items in \(A \bigcap B \bigcap C\) back into the final set.

Visually representing A, B, and C with a Venn diagram,

Project Euler Problem 1 Venn Diagram

To get back to the problem at hand, we can compute the size of these sets using the floor function. For example, the cardinality of A is:

$$ \mbox{card}(A) = \mbox{floor}\left( \frac{2001}{3} \right) = 667 $$
$$ \mbox{card}(B) = \mbox{floor}\left( \frac{2001}{4} \right) = 500 $$

Next, we subtract the duplicates (numbers with both A and B as factors):

$$ \mbox{card}(A \bigcap B) = \mbox{floor}\left( \frac{2001}{3 \cdot 4} \right) = 166 $$

Now subtract integers that have both a and c as multiples, or b and c as multiples:

$$ \mbox{card}(A \bigcap C) = \mbox{floor}\left( \frac{2001}{3 \cdot 5} \right) = 133 $$
$$ \mbox{card}(B \bigcap C) = \mbox{floor}\left( \frac{2001}{4 \cdot 5} \right) = 100 $$

And last but not least, those numbers with a, b, and c as factors were just removed twice, so we add them back in once:

$$ \mbox{card}(A \bigcap B \bigcap C) = \mbox{floor}\left( \frac{2001}{3 \cdot 4 \cdot 5} \right) = 33 $$

This gives a total number of multiples M below N with a or b as a factor, but not c:

$$ M = \mbox{floor}\left( \frac{N}{a} \right) + \mbox{floor}\left( \frac{N}{b} \right) - \mbox{floor}\left( \frac{N}{ab} \right) - \mbox{floor}\left( \frac{N}{ac} \right) - \mbox{floor}\left( \frac{N}{bc} \right) + \mbox{floor}\left( \frac{N}{abc} \right) $$

in our specific case,

$$ \begin{align} M &=& 667 + 500 - 166 - 133 - 100 + 33 \\ M &=& 801 \end{align} $$

Tags:    computer science    mathematics    factors    sequences    project euler   

Shortest Lattice Paths and Multiset Permutations

Posted in Mathematics

permalink

Table of Contents

The Lattice Paths Problem

I first came across the lattice paths problem in Project Euler problem 15. The question described a 2x2 square lattice, and illustrated the 6 ways of navigating from the top left corner to the bottom right corner by taking the minimum number of steps - 2 right steps and 2 down steps.

The question then asks for the number of minimum paths on a 20 x 20 grid. Needless to say, even without seeing the number, it should be obvious that enumerating all of these paths would get extremely expensive with grid dimensions growing beyond 10. That means this should be approached as an analytical combinatorics problem, not a computational combinatorics problem. As it turns out, there is a closed-form solution, and this is one of the few Project Euler questions that can be answered with the straightforward use of Wolfram Alpha (I solved it while boarding a bus).

But this is an interesting problem that goes beyond the Project Euler question - it has to do with a combinatorics problem that is maddeningly simple, yet surprisingly difficult to formulate - the problem of multisets.

Problem Formulation: Permutations

Thinking through the 2x2 grid, we already stated that the shortest path must consist of 2 right moves and 2 down moves - the number of paths is simply a question of the order in which these moves are made. Let us represent the path that moves right twice, then down twice, using the string

RRDD

Now we have a way of representing paths through the lattice - and we've turned our very specific lattice problem into a much more general combinatorics problem. How many unique permutations of the above path/string can we find?

Permutations of Unique Items

Let's suppose we have a string consisting of unique characters:

ABCDEFG

Now how many permutations of this string are there? The first letter can be any of the 7 characters, so we have 7 possibilities. The second letter can be any of the remaining 6 characters, so we have 7 * 6 possibilities. And so on down the line, until we get a total number of possible permutations of the string equal to

$$ 7! = 7 \times 6 \times 5 \times 4 \times 3 \times 2 \times 1 = 5040 $$

There are 5040 unique permutations of the string.

Our situation is complicated by the fact that some of our permutations are repeated. For example, if we label the two down moves as D1 and D2, we can choose the first move as D1 and the second move as D2, or the first move as D2 and the second move as D1 - the two are equivalent. This will eliminate some of the permutations.

Permutations of Items with Duplicates

Our case is slightly different: we have items with duplicates. This fits into a general combinatorics framework called stars-and-bars (link to Wikipedia article). In this framework, we are trying to determine the number of ways that we can partition a set of n objects into a set of \(k\) bins. The \(n\) objects are often denoted as star characters, and the k bins are formed by k-1 bar characters.

For example, if we are adding 5 components to a circuit board, and they can be any one of 9 possible components, we can represent this as the partitioning of 5 stars among 9 bins, or 8 bars. Here are some possibilities:

||||||||        <-- No choices made yet (9 bins, 8 partitions)
*|*||*||*|*||   <-- Mix of different components
||*|*||*|*||*   <-- Mix of different components
****|||*|||||   <-- Heavy on component 1
*|**|||**||||   <-- Two pairs

In the case of our lattice path problem, we have only two unique characters, D and R. We can think of the problem of generating permutations as inserting the down moves into a sequence of right moves. We have a certain number of locations where we can insert the down moves, and down moves can be inserted in order.

This is what's often called a stars-and-bars problem in combinatorics: trying to determine the number of permutations of items from multisets can be described as partitioning star characters with bar characters.

To formulate our problem in these terms, we can think of "distributing" our down moves as we move right through the lattice. On the 20x20 grid, we are going to make 20 right moves, and we can distribute our 20 down moves at any of 21 possible locations (columns of points on the lattice).

Thus, our n items, the 20 down moves, are being placed between the 20 right moves, which are the 20 bars that create 21 bins (21 locations to place the down moves).

Example

Considering the smaller 2x2 example, and replacing bars with R, we have, for a 2x2 lattice, six possibilities:

RR    <-- No choices made yet
**RR  <-- all on left
*R*R  <-- distributed... etc...
*RR*
R**R
R*R*
RR**

To enumerate the number of possibilities, we use the multichoose function, denoted "n multichoose k". This counts the number of ways to place n objects (D) into \(k\) bins (created by \(k-1\) other objects, R). The multichoose function is defined as (see https://en.wikipedia.org/wiki/Multiset#Counting_multisets wikipedia page for proper Latex notation - it looks like the binomial coefficient but with 2 parentheses):

$$ n \mbox{ multichoose } k = \binom{n+k-1}{n} $$

where, again, n is number of objects, being split into \(k\) partitions by \(k-1\) other objects.

Now, let's plug in the numbers for the 2 by 2 lattice. We get:

$$ n = 2, k = 3 $$

We are partitioning 2 down moves among 3 possible columns in the lattice. This gives:

$$ 2 \mbox{ multichoose } 3 = \binom{2+3-1}{2} = \binom{4}{2} = 10 $$

which is indeed the number of paths through the lattice.

Solving 2D Rectangular Lattice

To generalize, on a lattice of width \(W\) and height \(H\), we have \(W\) right moves that form \(W+1\) partitions, in which we are placing H items. The number of possible paths through the lattice is therefore equivalent to permutations of the string:

RRRR...(W times)...RRRDDDD...(H times)...DDDD

Now n and k are given by:

$$ n = H, k = W+1 $$

so the total number of possible paths through the W x H square lattice is:

$$ (H) \mbox{ multichoose } (W+1) = \binom{W+1+H-1}{H} = \binom{W+H}{H} $$

More Examples

The number of minimal paths through a 4 x 2 lattice (identical to the number of paths through a 2 x 4 lattice) is:

$$ P = \binom{4+2}{2} = \binom{4+2}{4} = 15 $$

The number of minimal paths through an 8x8 lattice is given by:

$$ P = \binom{8+8}{8} = 12,870 $$

and finally, the number of minimal paths through a 20 x 20 lattice is given by:

$$ P = \binom{20+20}{20} = 137,846,528,\mbox{XXX} $$

(This is the Project Euler problem 15 answer so last few digits are omitted.)

Solving 2D Square Lattice (Special Case)

If we use the above formulas for the special case where the dimensions of the grid are equal, such as the 20 x 20 case, we get the simpler and more symmetric formula:

$$ P = \binom{2D}{D} $$

where \(D\) is the dimension of the square grid.

Solving 3D Cuboid Lattice

(Note: cuboids are the 3D analogue of 2D rectangles.)

If we take this idea a step further, and use a slightly different combinatoric formula, we can generalize the problem to paths through higher dimensional lattices. This is an approach I came up with through trial and error, and some experimentation with Mathematica.

Suppose we have a 3D lattice, composed of 8 cubes, 2 on each side. Now we wish to know: how many shortest Manhattan distance paths are there through the lattice, from one corner to the other?

This can be re-cast, as we did above, as the problem of counting unique permutations of a string of the form:

UURRBB

where U denotes an up move, R denotes a right move, and B denotes a back move.

While we could use the multiset approach from above to describe this problem, it turns out that this approach is not powerful enough to describe the problem of arbitrary 3D lattices.

Let's pose the problem slightly more generally: we have C bags of moves or characters, each of a different size. We must use each character i precisely as many times as we have instances in its bag C_i. How many unique permutations are there?

Consider the following example string of length \(N = 14\), consisting of:

$$ \begin{array} N_x &=& 4 \mbox{ UP moves} \\ N_x &=& 4 \mbox{ RIGHT moves} \\ N_x &=& 6 \mbox{ BACK moves} \end{array} $$

These form a path through a 3D lattice of size \(4 \times 4 \times 6\):

UUUURRRRBBBBBB

The number of unique permutations can be computed by breaking this into sub-problems. Start by asking how many permutations there are of the string:

UUUUXXXXXXXXXX

(Treating the Rs and Bs as identical). Then we get:

$$ \binom{N}{N_x} = \binom{N_x + N_y + N_z}{N_x} = \binom{14}{10} = 1001 $$

Next, we deal with the other sub-problem, the Xs, by asking how many ways we can permute the following string:

RRRRBBBBBB

which is solved via another binomial coefficient. This number of permutations is given by:

$$ \binom{N_y + N_z}{N_y} = \binom{N_y + N_z}{N_z} = \binom{10}{4} = \binom{10}{6} = 210 $$

Now combining these, we get the overall number of permutations:

$$ P = \binom{14}{4} \cdot \binom{10}{4} = 210,210 $$

Solving 3D Cubic Lattice (Special Case)

If we have the special case of a perfect cubic lattice, the formula above reduces to the nice and symmetric:

$$ \dfrac{(3n)!}{(n!)^3} $$

Solving N-Dimensional Square Lattice (N-Dimensional Multisets)

Let's look at the example of a traversal of a 4D lattice, which we can think of as the evolution of a 3D traversal in time (a step in the fourth dimension would represent a "pause" in the 3D traversal).

Consider the traversal of a cube with dimensions \(3 \times 4 \times 5 \times 3\). Then

$$ N = 3 + 4 + 5 + 3 $$
$$ N_x = 3 $$
$$ N_y = 4 $$
$$ N_z = 5 $$
$$ N_t = 3 $$

A path on this 4D lattice has the form:

UUURRRRBBBBBWWW

(where \(W\) denotes wait, or a step in the time dimension).

The number of permutations is given by:

$$ P = \binom{N_x + N_y + N_z + N_t}{N_x} \cdot \binom{N_y + N_z + N_t}{N_y} \cdot \binom{N_z + N_t}{N_z} $$

For our specific example,

$$ P = \binom{3+4+5+3}{3} \cdot \binom{4+5+3}{4} \cdot \binom{5+3}{5} = 455 \cdot 495 \cdot 56 = 12,612,600 $$

Confirmed by Mathematica:

Permutations with Mathematica

Computing Square Roots: Part 2: Using Continued Fractions

Posted in Mathematics

permalink

Table of Contents

Continued Fractions

Let's start part 2 of our discussion of computing square roots by talking about continued fractions. When we first learn mathematics, we learn to count in the base 10 system: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9. We can construct representations of all of the integers using these 10 digits, by arranging them in a different order. So, for example, saying 125 is equivalent to saying:

$$ 125 = 1 \times 10^2 + 2 \times 10^1 + 5 \times 10^0 $$

Later on in our mathematical lives, we learn that we can use other integers as our base, or radix, by decomposing 125 into powers of that integer. For example, 125 can be decomposed into powers of 2. In a base 2 system, we have only two symbols, 0 and 1, so 125 can be represented as 1111101, which is equivalent to saying:

$$ 125 = 1 \times 2^6 + 1 \times 2^5 + 1 \times 2^4 + 1 \times 2^3 + 1 \times 2^2 + 1 \times 2 $$

Note that 1111101 is close to 1111111, which is \(2^7 = 128\).

As Knuth points out in his Art of Computer Programming, Part II:

"If our ancestors had invented arithmetic by counting with their two fists or their eight fingers, instead of their ten "digits," we would never have to worry about writing binary-decimal conversion routines. (And we would perhaps never have learned as much about number systems.)

This idea of an alternative radix to the traditional base 10 leads to entirely new number systems with their own interesting properties. However, it goes even further - as a high school student in 1955, Donald Knuth invented a number system with an imaginary, 4-symbol radix in base \(2i\) (where \(i = \sqrt{-1}\)), called the quater-imaginary base.

There are also number systems that use an irrational radix, such as phinary, which uses the golden ratio as its radix. Thus, the Golden Ratio \(\phi = \dfrac{1 + \sqrt{5}}{2}\) becomes \(\phi = 1\), \(2 = \phi^1 + \phi^-2\) becomes \(2 = 10.01\), \(5 = \phi^3 + \phi^{-1} + \phi^{-4}\) becomes \(5 = 1000.1001\), and so on. In the case of the irrational number \(\sqrt{5}\), this number becomes a rational number in base \(\phi\): \(\sqrt{5} = \phi^2 + \phi^{-1}\) becomes \(\sqrt{5} = 10.1\).

The single central idea behind these number systems is that the abstract set of integers \(\mathbb{Z}\) are being represented on an algebraic ring, which is a fundamental mathematical object. Initially the idea of a ring may seem strange, but it creates the foundations of modern number theory.

Continued Fraction Representations

(Note that while the next two sections will have a lot of "magic numbers," we will step through the procedure in the following sections, and it will be more clear where these numbers come from. There is also an excellent continued fractions calculator available online here.)

We can use other kinds of mathematical objects to represent different numbers. Another technique is to use continued fractions to represent numbers. These are well-studied mathematical objects that date back to Euclid's Elements. The basic idea is to create a recursive expression that involves fractions nested beneath other fractions:

$$ a_0+\cfrac{1}{a_1 +\cfrac{1}{a_2 +\cfrac{1}{ \begin{array}{@{}c@{}c@{}c@{}} a_3 + {}\\ &\ddots\\ &&{}+ \cfrac{1}{a_n} \end{array} }}} $$

For ease of writing, this is often written as:

$$ a_0 + \dfrac{1}{a_1 +} \quad \dfrac{1}{a_2 +} \quad \dfrac{1}{a_3 + } \quad \dots \quad \dfrac{1}{a_{n-1} + } \quad \dfrac{1}{a_n} $$

or, even shorter,

$$ [a_0; a_1, a_2, a_3, \dots] $$

These variables \(a_i\) are called the terms of the continued fraction. Continued fractions can be used to represent rational numbers, in which case the continued fraction representation terminates after a finite number of terms. For example, \(\dfrac{125}{3} = [41; 1, 2]\),

$$ \dfrac{125}{3} = 41 +\cfrac{1}{1 +\cfrac{1}{2}} $$

Continued fractions can also be used to represent irrational numbers, in which case the continued fraction representation is a repeating pattern of variable length. For example, \(\sqrt{14} = [3; \overline{1,2,1,6}]\), where the line over the last digits indicates that the pattern repeats infinitely as \(1, 2, 1, 6, 1, 2, 1, 6, 1, 2, 1, 6, \dots\):

$$ 3 + \cfrac{1}{ 1 + \cfrac{1}{ 2+\cfrac{1}{ 1+\cfrac{1}{ 6+\cfrac{1}{ \begin{array}{@{}c@{}c@{}c@{}} 1 + {}\\ &\ddots\\ &&{} \end{array} } } } } } $$

A few useful properties of these patterns, for square roots, are:

  • First, the integer portion (3 in the case above) is the largest integer whose square is less than the number (14).

  • Second, the sequence of integers that repeats is always palindromic, and it always begins repeating once it reaches a value of \(2 a_0\).

Here are a few more square roots represented as continued fractions, to help illustrate the above properties:

$$ \begin{array} _ \sqrt{19} &=& [4; \overline{2, 1, 3, 1, 2, 8}] \\ \sqrt{115} &=& [10; \overline{1, 2, 1, 1, 1, 1, 1, 2, 1, 20}] \\ \sqrt{988} &=& [31; \overline{2, 3, 4, 1, 20, 6, 1, 14, 1, 6, 20, 1, 4, 3, 2, 62}] \end{array} $$

Convergents

We will show an example of how to compute a continued fraction in the next section, but we will cover one additional topic first. In order to be useful, we need some way to evaluate the continued fraction representation. One way to do this is to repeatedly compute common denominators, perform the fraction addition, and inert the result. The fraction that results is equivalent to the continued fraction expression, but is a rational number and therefore easier to evaluate.

It turns out that this rational number is a very important quantity called the convergent of \(\sqrt{n}\). However, it is cumbersome to perform the operation of de-rationalizing the continued fraction. There is a useful shortcut that takes the form of a recurrence relation.

The \(n^{th}\) convergent is defined as the rational number that results when the continued fraction representation is carried out to \(n\) terms, then expanded and simplified into a rational number. If a number like \(\sqrt{19}\) has an infinite continued fraction representation, it means we can compute progressively more accurate rational approximations.

For example, we know from the above that \(\sqrt{19} = [4; \overline{2,1,3,1,2,8}]\). Using this fact, we can use successive terms to write successive linear approximations:

$$ \begin{array} _ \sqrt{19} &\approx& 4 \\ \quad &\approx& 4 + \frac{1}{2+0} \approx \frac{9}{2} \\ \quad &\approx& 4 + \cfrac{1}{2+\cfrac{1}{1}} \approx \frac{13}{3} \\ \quad &\approx& 4 + \cfrac{1}{2+\cfrac{1}{\cfrac{1}{3}}} \approx \frac{48}{11} \\ \end{array} $$

If we continue this sequence, we get a slew of approximations (there are also additional approximations between each of these terms...):

$$ \begin{array} _ \sqrt{19} &\approx& \frac{61}{14} \\ \quad &\approx& \frac{170}{39} \\ \quad &\approx& \frac{1421}{326} \\ \quad &\approx& \frac{3012}{691} \\ \quad &\approx& \frac{4433}{1017} \\ \sqrt{19} &\approx& \frac{16311}{3742} \end{array} $$

The numerator and denominator of the \(k^{th}\) convergent are denoted \(P_k\) and \(Q_k\), respectively, and can be computed through the recurrence relation:

$$ \dfrac{P_k}{Q_k} = \dfrac{a_k P_{k-1} + P_{k-2}}{a_k Q_{k-1} + Q_{k-2}} $$

where the initial values are \(P_{-1} = 1, P_{-2} = 0, Q_{-1} = 0, Q_{-2} = 1\) (making the first convergent equal to \(\frac{a_0}{1}\)). This can be used to compute successive approximations. Note that on a computational platform, you will quickly reach the end of your accuracy limit, so care must be taken for continued fraction sequences of longer than about 10 terms.

Example: Continued Fraction Coefficients of \(\sqrt{14}\)

Let's walk through an example of how to compute the continued fraction representation \([a_0; a_1, a_2, \dots, a_n]\) for a square root, and how to compute the \(k^{th}\) convergent.

Begin with the square root of a given number. For variety, we will use \(\sqrt{14}\). We begin by computing the nearest integer to 14's square root. If we don't have a square root routine, we can try integers by squaring them, and find the last integer that, when squared, is less than 14. This is \(a_0 = 3\). To interpret, \(a_0\) is the largest integer portion of our square root, with the residual portion \(\sqrt{14}-3\) the portion that will be represented by the continued fraction.

Here's what we have:

$$ \sqrt{14} = 3 + (\sqrt{14} - 3) \\ $$

We can drop the residual for an initial rational approximation of \(\sqrt{14} \approx 13\), but we can do better by going another step.

We recognize that the residual, which we want write in the form \(\dfrac{1}{\mbox{something}}\), can be written as \(\dfrac{1}{\dfrac{1}{\sqrt{14}-3}}\):

$$ \sqrt{14} = 3 + \dfrac{1}{ \frac{1}{\sqrt{14}-3} } $$

Now examine the inverse residual term \(\dfrac{1}{\sqrt{14}-3}\) for step 1, which we will call \(r_1\). Repeat the procedure we performed above: find the nearest integer portion to this quantity. In this case,

$$ \dfrac{1}{r_1} = \dfrac{1}{\sqrt{14}-3} = 1.34833... $$

Now split this into its integer portion, which becomes \(a_1\), and its remaining fractional portion, .34833...:

$$ a_1 = floor( \dfrac{1}{r_1} ) $$

and the new residual \(r\) is written in terms of the old residual and the coefficent \(a_i\):

$$ r_2 = \dfrac{1}{r_1} - a_1 $$

Generalizing, we get an iterative procedure to determine the coefficients \(a_i\):

$$ a_i = floor( \dfrac{1}{r_i} ) $$

followed by:

$$ r_{i+1} = \dfrac{1}{r_i} - a_i $$

where the initial values \(a_0, r_0\) are computed as mentioned above, and the rest of the values in the sequence follow.

Continuing for \(\sqrt{14}\), we get:

$$ \begin{array} _ a_0 &=& 3 \\ r_1 &=& \dfrac{1}{\sqrt{14}-3} = 1.348331 \\ a_1 &=& 1 \\ r_2 &=& \dfrac{1}{.348331} = 2.870829 \\ a_2 &=& 2 \\ r_3 &=& \dfrac{1}{0.870829} = 1.1483311 \\ a_3 &=& 1 \\ r_4 &=& \dfrac{1}{0.1483311} = 6.741676 \\ a_4 &=& 6 \end{array} $$

It is at this point that we see \(2 a_0\) and know that our (palindromic) sequence will repeat. (When we evaluate the convergents, we will utilize the palindromic nature of this sequence.)

Collecting these terms gives us the expected result: \(\sqrt{14} = [3; 1, 2, 1, 6]\).

This gives us an algorithmic procedure for computing the continued fraction representation \([a_0; a_1, a_2, a_n]\) of a number - with the important caveat, as mentioned above, that some integers have sequences that are extremely long before they repeat, making it impossible to find the full continued fraction representation without arbitrary precision libraries.

That being said, if you are only interested in the first 10 or so terms of the continued fraction representation, here is a static Java method to compute them:

    /** Find the (shorter than 10) continued fraction sequence for sqrt(n). 
     * This returns a list of integers, [a0, a1, a2, ...]
     *
     * @params n Number to compute the square root of.
     */
    public static List<Integer> continuedFractionSqrt(int n) {
        if(isPerfectSquare(n)) {
            throw new IllegalArgumentException("Error: n cannot be a perfect square.");
        }
        int niters = 10; // handbrake

        int ai = 0;
        double val = 0;
        double remainder = 0;
        List<Integer> coeffs = new ArrayList<Integer>();

        // Fencepost
        remainder = 1.0/Math.sqrt(n);

        for(int i=0; i<niters; i++) {
            val = 1.0/remainder;
            ai = floor(val);
            remainder = val - ai;
            coeffs.add(ai);
            if(coeffs.get(i) == 2*coeffs.get(0)) {
                break;
            }
        }
        return coeffs;
    }

    /** Check if x is a perfect square. */
    public static boolean isPerfectSquare(int x) { 
        int s = (int)(Math.round(Math.sqrt(x)));
        return x==(s*s);
    }

    /** Find the floor of a double. */
    public static int floor(double j) {
        return (int)(Math.floor(j));
    }

Here is a short program that uses this routine to compute the continued fraction representation of \(\sqrt{14}\):

public class SquareRootCF {
    public static void main(String[] args) {
        System.out.println("sqrt(14) = "+Convergents.continuedFractionSqrt(14));
        System.out.println("sqrt(19) = "+Convergents.continuedFractionSqrt(19));
    }
}

and the corresponding output:

$ javac SquareRootCF.java && java SquareRootCF
sqrt(14) = [3, 1, 2, 1, 6]
sqrt(19) = [4, 2, 1, 3, 1, 2, 8]

Example: Convergents of \(\sqrt{14}\)

We now turn to the task of computing the convergents of the continued fraction, which will yield successive rational numbers that are progressively better approximations to \(\sqrt{14}\).

We start with the expression given above for the \(k^{th}\) convergent:

$$ \dfrac{P_k}{Q_k} = \dfrac{a_k P_{k-1} + P_{k-2}}{a_k Q_{k-1} + Q_{k-2}} $$

with initial values \(P_{-1} = 1, P_{-2} = 0, Q_{-1} = 0, Q_{-2} = 1\). This yields the first "real" convergent:

$$ \dfrac{P_1}{Q_1} = \dfrac{a_0 + 0}{0 + 1} = 4 $$

Successive approximations will use the values \(P_1\) and \(Q_1\) to compute the next convergents.

$$ \dfrac{P_2}{Q_2} = \dfrac{P_1 a_1 + P_0}{Q_1 a_1 + Q_0} = \frac{11}{3} $$

Continuing in this fashion gives:

$$ \begin{array} \quad \dfrac{P_2}{Q_2} &=& \dfrac{15}{4} \\ \dfrac{P_3}{Q_3} &=& \dfrac{11}{3} \end{array} $$

and so on. This recurrence relation is easy to code up. It starts with the continued fraction coefficients for the given square root, and computes successive values of P and Q. The number of terms computed is specified by the user. Once it reaches the end of the sequence of continued fraction coefficients, it can start at the beginning again (the sequence is palindromic).

Finally, it returns the values of \(P_k\) and \(Q_k\), and of the successive convergents:

$$ \begin{array} _ \sqrt{14} &\approx& 3 \\ &\approx& 4 \\ &\approx& 11/3 \\ &\approx& 15/4 \\ &\approx& 101/27 \\ &\approx& 116/31 \\ &\approx& 333/89 \\ &\approx& 449/120 \\ &\approx& 3027/809 \\ &\approx& 3476/929 \\ &\approx& 9979/2667 \end{array} $$

Here is a static method in Java that will compute the convergents of a square root:

    /** 
     * Compute the convergents (rational representation of
     * continued fraction).
     *
     * This uses the recurrence relation:
     *
     * P_n     a_n P_n-1  + P_n-2
     * ---- = -----------------
     *  Q_n    a_n Q_n-1  + Q_n-2
     */
    public static long[] convergents(int n, int nterms) {
        long[] convergents = new long[2];

        List<Integer> cfrepr = continuedFractionSqrt(n);

        // Initial values for convergent recurrence relation
        long Pnm2 = 0; // P_{n-2}
        long Pnm1 = 1;
        long Qnm2 = 1;
        long Qnm1 = 0;
        long P = 0;
        long Q = 0;

        // Term 0 is the constant value a0.
        int accessindex = 0;
        for(int i=0; i<=nterms; i++) { 
            int an = cfrepr.get(accessindex);

            P = an * Pnm1 + Pnm2;
            Q = an * Qnm1 + Qnm2;

            Pnm2 = Pnm1;
            Pnm1 = P;

            Qnm2 = Qnm1;
            Qnm1 = Q;

            if(accessindex+1>=cfrepr.size()) { 
                // Ensure we keep repeating the sequence
                // if the sequence has fewer terms than
                // the user asks for.
                // This allows us to get really good
                // approximations for numbers.
                // This only works because the sequence
                // is palindromic.
                accessindex = 1;
            } else {
                accessindex++;
            }
        }

        convergents[0] = P;
        convergents[1] = Q;

        return convergents;
    }

Here is a simple driver program that prints out several convergents for \(\sqrt{14}\):

public class SquareRootCF {
    public static void main(String[] args) {
        for(int i=1; i<=10; i++) {
            long[] conv = Convergents.convergents(14,i);
            System.out.println("Convergent "+i+": "+conv[0]+"/"+conv[1]);
        }
    }
}

and the corresponding console output:

$ javac SquareRootCF.java && java SquareRootCF
Convergent 1: 4/1
Convergent 2: 11/3
Convergent 3: 15/4
Convergent 4: 101/27
Convergent 5: 116/31
Convergent 6: 333/89
Convergent 7: 449/120
Convergent 8: 3027/809
Convergent 9: 3476/929
Convergent 10: 9979/2667

References

  1. Hardy, G. H. A Course of Pure Mathematics. Cambridge University Press, Tenth Edition (1908-1967).

  2. Knuth, Donald. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. Addison-Wesley Publishing Company, Second Edition (1975).