Continued Fractions: Introduction and Gosper Arithmetic
Continued Fractions (Khinchin)
Download the notebook
Continued fractions are expressions of the form \(\lbrack a_0;a_1,a_2,\ldots\rbrack = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{\cdots}}}.\) Every real number has a unique expression as a continued fraction, so long as we exclude representations which end in 1, since \(a_0;a_1,\ldots,a_n,1] = [a_0;a_1,\ldots,a_n+1].\) Note that this is no worse (and probably quite a bit better) than the ambiguity in decimal expansions, for which \(1 = 0.999\cdots.\)
There are a few reasons why continued fractions are much better behaved than decimal (or general base \(b\)) expansions. All of them boil down to the fact that a continued fraction representation of \(x\) is intrinsic to \(x\), while a base \(b\) expansion forces \(x\) to wear a stranger’s clothes.
Rational Numbers
In a decimal expansion, rational values can have infinite expansions: \(\begin{align} \frac{1}{3} &= 0.3333\cdots = \sum_{n=1}^\infty 3\cdot 10^-n,\\ \frac{1}{99} &= 0.12345\cdots = \sum_{n=1}^\infty n\cdot 10^{-n},\\ \frac{1}{89} &= 0.011235 = \sum_{n=1}^\infty F_{n-1} \cdot 10^{-n}, \end{align}\) where \(F_n\) are the Fibonacci numbers starting with \(F_0 = 0\), \(F_1 = 1\). This last fact is fun, but it has nothing to do with \(\frac{1}{89}\), except for the coincidence that \(89 = 10^2 - 10 - 1\). You can check that \(\frac{1}{239} = \frac{1}{16^2-16-1}\) also has Fibonacci numbers in its hexadecimal expansion.
On the contrary, a number has a finite continued fraction representation if and only if it is rational. Clearly no irrational number is a finite continued fraction; to establish the converse, let \(\frac{p}{q}\) be a rational number with \(p<q\) and \(\gcd(p,q) = 1\). Then the continued fraction expansion for \(\frac{p}{q}\) begins \(\frac{p}{q} = \frac{1}{\lfloor q/p\rfloor + \frac{1}{r}}.\) Solving for \(r\), we find \(r = \frac{1}{q/p - \lfloor q/p \rfloor} = \frac{p}{q-p\lfloor q/p\rfloor}.\) We have written \(r\) as a quotient of integers with denominator \(q-p\lfloor q/p\rfloor\), which (since \(p<q\)) is strictly less than \(q\). Reducing \(r\) to lowest terms can only lower the denominator further, so it follows that the denominator of \(r\) is less than \(q\). We may then apply the same procedure to \(r\), obtaining a new remainder \(r'\) with still smaller denominator, and so on; the process will at some point terminate when the denominator becomes 1, that is, when we reach an integer.
Quadratic Irrationals
In addition to representing rational numbers well, there is a remarkable pattern in the continued fractions of quadratic irrationals. As a simple example, let’s start with \(\sqrt{2}\). As a continued fraction, it is \(\sqrt{2} = [1;2,2,2,\cdots].\) It is simple to prove this. We know that \(\sqrt{2}\) is a solution to \(x^2 = 2\), which we can write as \(x^2 - 1 = 1\), or better yet, \(x = 1+\frac{1}{1+x}.\) Now let \(x = [1;a_1,a_2,\cdots]\). This equation then reads \(1;a_1,a_2,\cdots] = 1+\frac{1}{[2;a_1,a_2,\cdots]} = [1;2,a_1,a_2,\cdots].\) Comparing the continued fraction term by term, we find that \(a_i = 2\) for all \(i=1,2,\cdots\).
Let’s look at some other quadratic irrationals: \(\begin{align} \sqrt{3} &= [1;1,2,1,2,1,2,\cdots], \\ \sqrt{23} &= [4;1,3,1,8,1,3,1,8,\cdots],\\ \varphi = \frac{1+\sqrt{5}}{2} &= [1;1,1,1,\cdots]. \end{align}\) They’re all periodic! This is quite astonishing from a decimal perspective; there’s no way to look at the digits of an irrational number like \(\sqrt{2} = 1.4142135\ldots\) and recognize that it is the solution to a quadratic equation over rationals, but continued fractions make this property manifest.
The general fact here is that a continued fraction is eventually periodic if and only if it represents a quadratic irrational number. The proof is a bit technical, but the idea is to establish a recurrence among the terms of the continued fraction, the same way we did for \(\sqrt{2}\). It quickly gets tricker to construct the recurrence explicitly, though. For example, take \(\sqrt{3}\). We start with \(x^2 = 3\), and then write this as \(x^2 + x - 2 = x + 1\). This weird rearrangement allows us to factor the left hand side as \((x-1)(x+2)\), and then we express \(x+2 = (x+1)\left(1+\frac{1}{x+1}\right)\), so that we have \((x-1)(x+1)\left(1+\frac{1}{2+(x-1)}\right) = x+1.\) Canceling the factor of \(x+1\), and rearranging a bit more, we find \(x = 1 + \frac{1}{1 + \frac{1}{2 + (x-1)}},\) which gives the recurrence we’re looking for: \([1;a_1,a_2,\cdots] = [1;1,2,a_1,a_2,\cdots]\).
Rational Approximation
Let’s look at \(\pi\) in both decimal expansion and continued fraction expansion: \(\begin{align} \pi &= 3.14159235\cdots \\ &= [3;7,15,1,292,1,1,2,1,\cdots] \end{align}\) Both expansions are infinite, so to report the value of \(\pi\), we have to truncate either one. If we take one digit of either representation, we just get 3; what if we take two digits? In decimal, we have 3.1, but in continued fractions, we have \(3 + \frac{1}{7} = \frac{22}{7} = 3.1429\cdots,\) which is a familiar rational approximation to \(\pi\) and also about 40 times closer to the actual value than 3.1. How about three digits? Then we have 3.14 versus \(3 + \frac{1}{7 + \frac{1}{15}} = \frac{333}{106} = 3.1415094\cdots,\) and the continued fraction is about 20 times closer.
This is a general phenemonon: truncating a continued fraction generates a convergent of the number it represents. A convergent to \(\alpha\) is a rational number \(\frac{p}{q}\) such that, for any other \(\frac{p'}{q'}\) with \(q' \le q\), \(|q\alpha - p| \le |q'\alpha - p'|.\) Every convergent to \(\alpha\) is a truncation of its continued fraction representation, and every truncation of a continued fraction representation is a convergent.
Better yet, we can say just how close the convergents get to the value of \(\alpha\). Let \(\alpha = [a_0;a_1,a_2,\ldots]\), and let the \(n\)th convergent be \(\frac{p_n}{q_n}\), with \(p_0 = a_0\) and \(q_0 = 1\). Then \(\frac{1}{(a_{n+1}+2)q_n^2} < \left|\alpha-\frac{p_n}{q_n}\right| < \frac{1}{a_{n+1}q_n^2}.\) This formula shows that large numbers in a continued fraction expansion imply that the prior convergent is getting close to the actual value. For instance, \([3;7,15,1,292] = \frac{355}{113} = 3.14159292\cdots\) is quite close to \(\pi\), and in fact differs by no more than \(\frac{1}{292(113)^2}\). In this case the bound is rather tight; the convergent differs by \(2.67\times 10^{-7}\), and the bound is \(2.68\times 10^{-7}\).
On the contrary, remember that the golden ratio is \(\varphi = [1;1,1,1,\cdots]\). This means that \(\varphi\) is the most difficult number to approximate by rationals, since \(a_{n+1}q_n^2\) grows as slowly as possible. And in fact, you already know its convergents:
Convergents[(1+Sqrt[5])/2,10]
This pattern follows from a simple recurrence relation obeyed by the convergents: \(\begin{align} p_n &= a_n p_{n-1} + p_{n-2},\\ q_n &= a_n q_{n-1} + q_{n-2}. \end{align}\)
Transcendentality
It is tempting to think that there might be some tower of properties of the continued fractions of algebraic numbers: rational \(\leftrightarrow\) finite, quadratic \(\leftrightarrow\) periodic, cubic \(\leftrightarrow\) ? But unfortunately, little is known about the continued fractions for algebraic numbers of higher degree.
However, we can say something quite strong about the rational approximation of algebraic numbers. Roughly speaking, the higher the degree of an algebraic number, the better the rational approximations it admits. The precise statement is the converse:
Lemma: Let \(\alpha\) be an irrational number such that \(f(\alpha) = 0\), where \(f\) is a polynomial with integer coefficients and degree \(n\). Then there exists \(C > 0\) such that for any rational number \(\frac{p}{q}\), \(\left|\alpha-\frac{p}{q}\right| \ge \frac{A}{q^n}\)
For instance, let \(\alpha\) be a quadratic irrational. We have already shown that convergents come no closer than \(\frac{1}{(a_{n+1}+2)q_n^2}\), and for a quadratic irrational, the numbers \(a_n\) are periodic and thus bounded, so we can choose \(A = \frac{1}{\max a_n}\) to satisfy the requirement of the lemma.
Now, what if we construct a number such that, for every \(n > 0\), there are infinitely many rationals \(\frac{p}{q}\) for which \(\left|\alpha - \frac{p}{q}\right| < \frac{1}{q^n}.\) Then for any \(A\) and \(n\), we can pick an exponent \(m = n + r\) such that \(2^{-r} < A\), and then this condition violates the statement of the lemma. Thus, such a number is not algebraic of any degree \(n\), and so it must be transcendental.
How could we construct such a number? It’s relatively easy to do it with decimal expansions containing long strings of zeroes, but this lacks a certain panache. Continued fractions have all the pinache. All we have to do is write a continued fraction term by term, always taking \(a_n \ge q_{n-1}^n\). Then the \(n\)th convergent satisfies \(\left|\alpha - \frac{p}{q}\right| \le \frac{1}{a_{n+1} q_n^2} \le \frac{1}{q_n^{n+3}},\) and so \(\alpha\) satisfies the desired property, and is hence transcendental.
Numbers having this approximability property are called Liouville numbers, and this was the argument Liouville used to first establish the existence of transcendental numbers. Nowadays this can be proved in a single breath, since algebraic numbers are countable and real numbers are uncountable, but Liouville’s explicit construction is nice.
Arithmetic
Continued fractions are so beautiful – why don’t we use them more often? You might object that they’d make arithmetic more difficult; it’s not at all obvious how to add or multiply two continued fractions. In fact, it’s not even obvious how to double a continued fraction or add \(\frac{2}{3}\) to it.
But non-obvious does not mean impossible. Sometime in the 1970s, Bill Gosper devised an algorithm for computing term by term with continued fractions.
Gosper’s method takes some getting used to, but at the end of the day I’d argue it’s at least as intuitive than typical arithmetic-by-hand methods used on decimal expansions. To compute the decimal expansion of \(\alpha\), where let’s say \(1 \le \alpha < 10\) for simplicity, one is tasked with computing \(d_0 = \alpha\pmod{10}\), and then recursing on \(10(\alpha-d)\). To compute the continued fraction expansion of \(\alpha\), one instead computes \(a_0 = \lfloor \alpha \rfloor\) and recurses on \(\frac{1}{\alpha-a_0}\).
For instance, let’s compute \(\frac{17}{6}\) as a continued fraction. Its integer part is 2, so \(a_0 = 2\), and now we have to compute \(\frac{1}{17/6 - 2} = \frac{6}{5}\). Its integer part is 1, so now we have to compute \(\frac{1}{6/5-1} = 5\). But 5 is an integer, so we’re done: \(\frac{17}{6} = [2; 1, 5].\)
This is fine for expanding a known rational value into a continued fraction, but what about operating on continued fractions themselves? For instance, what if we want to compute \(4\varphi\) as a continued fraction? Note that this is a totally different game from computing \(4\varphi\) as a decimal. In a decimal expansion, we’d be forced to truncate and then multiply by 4 to get an approximation to \(4\varphi\). With continued fractions, we know \(\varphi\) exactly, and we can compute \(4\varphi\) exactly.
To do so, we follow the same process as before: find the integer part, and then recurse. To find the integer part, we’ll use the following fun fact: if \(a,b,c,d\ge 0\), then the mediant \(\frac{a+c}{b+d}\) lies between \(\frac{a}{c}\) and \(\frac{b}{d}\). So, let’s write \(4\varphi\) as \(\frac{4\varphi + 0}{0\varphi + 1}.\) This is the mediant of \(\frac{4}{0}\) and \(\frac{0}{1}\). Alas, there are many integers between 0 and \(\infty\), so we don’t know the integer part yet. But we know that the integer part of \(\varphi\) is 1, so we can write \(\varphi = 1 + \frac{1}{r_0}\), and rewrite our fraction as \(\frac{4\varphi + 0}{0\varphi + 1} = \frac{4(1+1/r_0) + 0}{0(1+1/r_0) + 1} = \frac{4r_0 + 4}{r_0 + 0}.\) Now we have the mediant of \(\frac{4}{1}\) and \(\frac{4}{0}\). Not quite what we need, but we know \(r_0 = 1 + \frac{1}{r_1}\), so this gives \(\frac{4r_0 + 4}{r_0 + 0} = \frac{4(1+1/r_1) + 4}{1 + 1/r_1 + 0} = \frac{8r_1+4}{r_1+1}.\) Now this is the mediant of \(\frac{8}{1}\) and \(\frac{4}{1}\), so the integer part must be either 4 or 5. One more step: \(\frac{8r_1+4}{r_1+1} = \frac{8(1+1/r_2) + 4}{1+1/r_2 + 1} = \frac{12r_2 + 8}{2r_2+1}.\) Now we know we’re between 6 and 8; so close! On the next step we get \(\frac{20r_3 + 12}{3r_3+2}\), and this shows that the integer part is 6. So we write down 6 and then work on computing \(\left(\frac{20r_3+12}{3r_3+2}-6\right)^{-1} = \frac{3r_3+2}{2r_3}.\) This expands in two steps to \(\frac{8r_5+5}{4r_5+2}\), which shows that the next digit is 2. Then we subtract 2, invert, and work on computing \(\frac{4r_5+2}{1}\). This eventually becomes \(\frac{26r_9+16}{3r_9+2}\), which gives an 8, and then we recurse and compute \(\frac{3r_9+2}{2r_9}\). But wait! Since \(r_9 = r_3\) (by known periodicity of \(\varphi\)), we’ve already computed this, and so we know that the \(2, 8\) repeats: \(4\varphi = [6; 2,8,2,8,\cdots].\)
To recap: this method allows us to compute term by term the continued fraction expansion for any expression \(\frac{ax+b}{cx+d}\) in terms of the expansion of \(x\). Better yet, if \(x\) is known exactly, either as a finite continued fraction or one with a repeating part, we can compute \(\frac{ax+b}{cx+d}\) exactly in the same sense. Better yet still, we compute the function of \(x\) term by term via accessing \(x\) term by term, so this algorithm can be “fed into itself” in a natural way, allowing us to chain computations together.
Let’s put this together into some actual code. We’ll focus on rational numbers first. We need a symbol to convert rational numbers to continued fractions, as in the example above with \(\frac{17}{6}\). This is effectively just Euclid’s algorithm for the GCD, where we store quotients as terms in the continued fraction.
CF::usage =
"CF[a,b,...] represents the continued fraction a + 1/(b + 1/(...))";
ToCF[x_Integer] := CF[x];
ToCF[Rational[p_, q_]] :=
If[q == 0, CF[], Join[CF[Floor[p/q]], ToCF[q/Mod[p, q]]]];
ToCF[17/6]
Gosper calls functions of the form \(\frac{ax+b}{cx+d}\) “homographic,” so that can be our Mathematica symbol for these functions. We can write our logic for computing homographic functions as follows:
Homographic[{{a_Integer, b_Integer}, {c_Integer, d_Integer}}, x_CF] :=
Which[
Length[x] == 0,(*x is CF[], which we can think of as infinity*)
ToCF[a/c],
a == c == 0, (*does not depend on x*)
ToCF[b/d],
c != 0 && d != 0 && Floor[a/c] == Floor[b/d], (*we know the integer part*)
Join[CF[Floor[a/c]], Homographic[{{c, d}, {Mod[a, c], Mod[b, d]}}, x]],
True, (*rewrite using a term of x*)
Homographic[{{a First[x] + b, a }, {c First[x] + d, c }}, Rest[x]]
];
For example, to double \(\frac{17}{6}\), we use
Homographic[{{2, 0}, {0,1}}, ToCF[17/6]]
ToCF[17/3]
Of course, this still doesn’t answer the original question: how to add, subtract, multiply, or divide two continued fractions. For this we have to generalize to what Gosper calls “bihomographic” functions, \(\frac{a xy + b x + c y + d}{e xy + f x + g y + h}.\) If \(\frac{a}{e}, \frac{b}{f}, \frac{c}{g}, \frac{d}{h}\) all have the same integer part, then we can write down this integer part and recurse, just as before. If they don’t all have the same integer part, then we have to think slightly harder, because we could expand using either a term of \(x\) or a term of \(y\). Just about any policy (even random choice) would get to the answer eventually, but one way to do it is to expand using \(x\) when \(\left|\frac{b}{f}-\frac{d}{h}\right| > \left|\frac{c}{g}-\frac{d}{h}\right|\). The full logic is as follows.
BiHomographic[{{a_Integer, b_Integer, c_Integer,
d_Integer}, {e_Integer, f_Integer, g_Integer, h_Integer}}, x_CF,
y_CF] :=
Which[
Length[x] == 0,(*x is infinity*)
Homographic[{{a,b},{e,f}},y],
Length[y] == 0,(*y is infinity*)
Homographic[{{a,c},{e,g}},x],
a == b == e == d == 0,(*does not depend on x*)
Homographic[{{c, d}, {g, h}}, y],
a == c == e == f == 0,(*does not depend on y*)
Homographic[{{b, d}, {f, h}}, x],
e != 0 && f != 0 && g != 0 && h != 0 &&
Floor[a/e] == Floor[b/f] == Floor[c/g] ==
Floor[d/h],(*we know the integer part*)
Join[CF[Floor[a/e]],
BiHomographic[{{e, f, g, h}, Mod[{a, b, c, d}, {e, f, g, h}]}, x,
y]],
f g h == 0 || Abs[b/f - d/h] > Abs[c/g - d/h],(*expand x*)
BiHomographic[{{a First[x] + c, b First[x] + d, a,
b }, {e First[x] + g, f First[x] + h, e, f}}, Rest[x], y],
True,(*expand y*)
BiHomographic[{{a First[y] + b, a, c First[y] + d,
c }, {e First[y] + f, e, g First[y] + h, g}}, x, Rest[y]]
];
This one algorithm allows us to add, subtract, multiply, and divide continued fractions. We can tell Mathematica this fact explicitly:
Plus[x_CF, y_CF] ^:= BiHomographic[{{0, 1, 1, 0}, {0, 0, 0, 1}}, x, y];
Times[x_CF, y_CF] ^:=
BiHomographic[{{1, 0, 0, 0}, {0, 0, 0, 1}}, x, y];
CF /: Times[x_CF, -1] := Homographic[{{-1, 0}, {0, 1}}, x];
CF /: Power[x_CF, -1] := Homographic[{{0, 1}, {1, 0}}, x];
Now we can do arithmetic directly with the CF expressions.
ToCF[1/2]+ToCF[1/5]
ToCF[1/2+1/5]
ToCF[34/3] ToCF[56/5]
ToCF[(34/3)(56/5)]
A few more definitions allows us to mix in regular rational numbers:
CF /: Plus[x_CF, y_Integer] := Homographic[{{1, y}, {0, 1}}, x];
CF /: Times[x_CF, y_Integer] := Homographic[{{y, 0}, {0, 1}}, x];
CF /: Plus[x_CF, Rational[p_, q_]] := Homographic[{{q, p}, {0, q}}, x];
CF /: Times[x_CF, Rational[p_, q_]] := Homographic[{{p, 0}, {0, q}}, x];
So, for instance:
eighth = ToCF[1/8];
Table[{i eighth, ToCF[i/8]}, {i, 8}]//TableForm
Finally, for convenience, let’s allow Normal to unwrap a continued fraction:
Normal[CF[a_, b__]] ^:= a + 1/Normal[CF[b]];
Normal[CF[a_]] ^:= a;
For example, we can use this to see what happens if we square \([1;2,\cdots,2]\), a convergent of \(\sqrt{2}\):
Normal/@Table[CF[1,Sequence@@ConstantArray[2,i]] CF[1,Sequence@@ConstantArray[2,i]] - 2,{i,10}]
Gosper Arithmetic with Periodic Fractions
It would be nice to do this calculation exactly and get 2. We saw how it’s possible to compute \(4\varphi\) exactly as a continued fraction, by recognizing when the algorithm begins to repeat itself. To write a continued fraction with a repeating part, we’ll put the repeating part at the end as a list, and tell Mathematica how to interpret this:
Rest[CF[{a__}]] ^:= CF[RotateLeft[{a}]]
First[CF[{a__}]] ^:= First[{a}];
We can tell Mathematica to check when the algorithm repeats itself by passing a cache variable through the recursion. First, with the homographic functions:
ClearAll[Homographic];
Homographic[{{a_Integer, b_Integer}, {c_Integer, d_Integer}}, x_CF] :=
Homographic[{{a, b}, {c, d}}, x, <||>];
Homographic[{{a_Integer, b_Integer}, {c_Integer, d_Integer}}, x_CF,
cache_] := Which[
KeyExistsQ[cache, {{a, b}, {c, d}, x}],
CF[cache[{{a, b}, {c, d}, x}]],
Length[x] == 0,
ToCF[a/c],
a == 0 && c == 0,
ToCF[b/d],
c != 0 && d != 0 && Floor[a/c] == Floor[b/d],
Join[CF[Floor[a/c]],
Homographic[{{c, d}, {Mod[a, c], Mod[b, d]}}, x,
Append[ Join[#, {Floor[a/c]}] & /@
cache, {{a, b}, {c, d}, x} -> {Floor[a/c]}]]],
True,
Homographic[{{a First[x] + b, a }, {c First[x] + d, c}}, Rest[x],
cache]
];
For instance, we can repeat the calculation of \(4\varphi\) above:
phi = CF[{1}];
4 phi
This is technically correct, but let’s tell Mathematica how to simplify it:
Simplify[CF[a___, b__, {b__}]] ^:= CF[a, {b}];
Simplify[CF[a___, {b__, b__}]] ^:= CF[a, {b}];
4 phi // Simplify
Much better. We can also try a more complicated example, and use Mathematica’s built-in ContinuedFraction symbol to check that things are working properly:
Simplify[4/3+CF[1,{2}]]
ContinuedFraction[4/3+Sqrt[2]]
Excellent. Now what about bihomographic functions, for instance to show that \([1;2,2,\cdots]^2 = 2\)? We could add a similar cache to the BiHomographic symbol, but there are two issues. First, to compute \([1,2,2,\cdots]^2\), we would first need its integer part – but since the answer is an integer, this amounts to the full computation. Thus, for any finite number of terms we extract from the two factors, we’ll never be able to nail down that the product is 2; we’ll always have the mediant of fractions slightly above and slightly below 2. Second, even in cases where the algorithm can generate terms of the continued fraction (such as for multiplying \(\sqrt{2}\) and \(\sqrt{3}\), getting \(\sqrt{6} = [2;2,4,2,4,\cdots]\)), the fractions appearing in Gosper’s algorithm are not guaranteed to repeat, and so the algorithm may not be able to recognize a pattern.
I suspect both of these issues can be solved by looking at recurrence equations on the coefficients \(a,\ldots,h\), but this is a problem for another time.
Enjoy Reading This Article?
Here are some more articles you might like to read next: