Continued fractions in Python

August 16, 2023 at 10:40 AM by Dr. Drang

My last post ends with “one last thing” about continued fractions. That turned out to be a lie. After playing around a bit more, I decided I should have some functions that compute continued fractions in Python, so I looked around for continued fraction libraries. I found some, but none of them seemed like the library. So I decided to write my own.

Let’s review some notation and properties. A continued fraction is one in which the denominator contains a fraction, and the denominator of that fraction contains a fraction, and so on.

x = a_{0} + \frac{1}{a_{1} + \frac{1}{a_{2} + \frac{1}{a_{3} + \dots}}}

This is considered the standard form for continued fractions, where the numerators are all ones. You can write out a continued fraction with other numbers as the numerators, but it can always be reduced to this form.

If $x$ is a rational number, then the continued fraction has a finite number of terms and will end with a $1 / a_{n}$ term. If $n = 4$ , for example, the fraction will look like this:

x = a_{0} + \frac{1}{a_{1} + \frac{1}{a_{2} + \frac{1}{a_{3} + \frac{1}{a_{4}}}}}

If $x$ is irrational, the continued fraction has an infinite number of terms, although terms may repeat. Famously, the golden ratio goes on forever and all the terms are one:

ϕ = 1 + \frac{1}{1 + \frac{1}{1 + \frac{1}{1 + \dots}}}

A less explicit but far more compact way to display a continued fraction is to just show the $a$ terms as a bracketed list:

x = [a_{0}; a_{1}, a_{2}, a_{3}, \dots]

It’s common to use a semicolon to separate the $a_{0}$ term from the others. Mathematica doesn’t do that because it’s more convenient to just use a list. As we saw in the last post, the first five terms of the continued fraction for $π$ is

In[1]:= ContinuedFraction[Pi, 5]

Out[1]= {3, 7, 15, 1, 292}

where Mathematica uses braces to surround its lists. We’ll use this same idea in Python, where the lists are bracketed.

A segment of a continued fraction, $s_{k}$ , is a finite continued fraction consisting of the first $k + 1$ terms of $x$ :

s_{k} = [a_{0}; a_{1}, a_{2}, \dots, a_{k}]

A remainder, $r_{k}$ , is all the terms starting with the $k^{t h}$ and continuing on, whether the continued fraction is finite or infinite:

r_{k} = [a_{k}; a_{k + 1}, a_{k + 2}, \dots]

So any continued fraction can be broken into a segment, $s_{k - 1}$ , and a remainder, $r_{k}$ .

A convergent is the rational number corresponding to a segment. Convergents are what we use to get rational approximations of numbers. In the last post, we did this

In[2]:= Convergents[ContinuedFraction[Pi, 5]]

            22  333  355  103993
Out[2]= {3, --, ---, ---, ------}
            7   106  113  33102

to see why $22 / 7$ and $355 / 113$ are good rational approximations of $π$ .

An interesting property of convergents is that those from even-indexed segments— $s_{0}$ , $s_{2}$ , and so on—bound $x$ from below, and those from odd-indexed segments— $s_{1}$ , $s_{3}$ , and so on—bound $x$ from above. (If $x$ is rational, then there is a final convergent which is equal to $x$ , regardless of whether it’s even or odd.) The even convergents form an increasing sequence; the odd convergents form a decreasing sequence.

OK, if you want more you can go to the Wikipedia page or get a copy of Khinchin’s book¹. Let’s move onto the code.

The function I wrote, continued, returns a tuple of

the continued fraction of the argument as a list of integers; and
the convergents of the argument as a list of Fractions.

Fractions are a Python type supplied by the fractions library.

continued is the only function in cfractions.py, a file I’ve saved in my site-packages directory. This makes it easy to import when I’m working in Jupyter:

In [1]: import math

In [2]: from cfractions import continued

In [3]: continued(math.pi)
Out[3]:
([3, 7, 15, 1, 292],
 [Fraction(3, 1),
  Fraction(22, 7),
  Fraction(333, 106),
  Fraction(355, 113),
  Fraction(103993, 33102)])

Here’s the code:

python:
 1:  from fractions import Fraction
 2:  from math import isclose
 3:  
 4:  def continued(x, terms=20, rel_tol=1e-9, abs_tol=0.0):
 5:    'Return the continued fraction and convergents of the argument.'
 6:    # Initialize, using Khinchin's notation
 7:    a = []       # continued fraction terms
 8:    p = [0, 1]   # convergent numerator terms (-2 and -1 indices)
 9:    q = [1, 0]   # convergent denominator terms (-2 and -1 indices)
10:    s = []       # convergent terms
11:    remainder = x
12:  
13:    # Collect the continued fraction and convergent terms
14:    for i in range(terms):
15:      # Compute the next terms
16:      whole, frac = divmod(remainder, 1)
17:      an = int(whole)
18:      pn = an*p[-1] + p[-2]
19:      qn = an*q[-1] + q[-2]
20:      sn = Fraction(pn, qn)
21:  
22:      # Add terms to lists
23:      a.append(an)
24:      p.append(pn)
25:      q.append(qn)
26:      s.append(Fraction(sn))
27:  
28:      # Convergence check
29:      if isclose(x, float(sn), rel_tol=rel_tol, abs_tol=abs_tol):
30:        break
31:  
32:      # Get ready for next iteration
33:      remainder = 1/frac
34:  
35:    # Return the tuple of the continued fraction and the convergents
36:    return(a, s)

The terms of the continued fraction are calculated using a form of Euclid’s algorithm for finding the greatest common divisor (GCD) of two numbers. The numerators and denominators of the convergents are calculated using the recurrence relations,

p_{k} = a_{k} p_{k - 1} + p_{k - 2}

q_{k} = a_{k} q_{k - 1} + q_{k - 2}

You may be wondering why I’m calculating both the continued fraction and the convergents in the same function instead of doing them separately as Mathematic does. Two reasons:

First, I typically want both the continued fraction and its convergents, so there’s no point in forcing me to make two function calls in the usual case.
Second, I want to be able to set specific convergence criteria, and to do that I need to know how close I am to the input argument as new terms of the continued fraction are calculated. That means I need to calculate the convergents along with the continued fraction terms.

The three optional parameters to the function, terms, rel_tol, and abs_tol, set the convergence criteria. terms is an upper bound on the number of continued fraction terms that will be calculated, no matter what the other tolerance values are. rel_tol and abs_tol, are relative and absolute tolerance values that can stop the process before the terms limit is reached. Their names and default values are taken from the isclose function of the math library, which is used on Line 29. For example, we could set an absolute tolerance on our rational estimate of $π$ this way:

In [4]: continued(math.pi, rel_tol=0, abs_tol=1e-12)
Out[4]:
([3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3],
 [Fraction(3, 1),
  Fraction(22, 7),
  Fraction(333, 106),
  Fraction(355, 113),
  Fraction(103993, 33102),
  Fraction(104348, 33215),
  Fraction(208341, 66317),
  Fraction(312689, 99532),
  Fraction(833719, 265381),
  Fraction(1146408, 364913),
  Fraction(4272943, 1360120)])

We’ve hit our tolerance because

π - \frac{4272943}{1360120} = 4 \times 10^{- 13}

I like the idea of having control over the convergence criteria. Mathematica’s second argument to ContinuedFraction gives the equivalent of my terms parameter, but its precision control is, as far as I can tell, entirely internal—there’s no way for the user to set a tolerance.

On the other hand, a disadvantage of my function is that its precision is limited to that of Python floats, whereas Mathematica will give you as much precision as you as for. I can’t, for example, ask for an abs_tol of $1 \times 10^{- 20}$ and expect to get a correct answer:

In [5]: continued(math.pi, rel_tol=0, abs_tol=1e-20)
Out[5]:
([3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3],
 [Fraction(3, 1),
  Fraction(22, 7),
  Fraction(333, 106),
  Fraction(355, 113),
  Fraction(103993, 33102),
  Fraction(104348, 33215),
  Fraction(208341, 66317),
  Fraction(312689, 99532),
  Fraction(833719, 265381),
  Fraction(1146408, 364913),
  Fraction(4272943, 1360120),
  Fraction(5419351, 1725033),
  Fraction(80143857, 25510582),
  Fraction(245850922, 78256779)])

Mathematica will happily use as many digits as needed and do so correctly. It tells us that Python screwed up in the 14th term:

In[3]:= ContinuedFraction[Pi, 15]

Out[3]= {3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1}

I’m not especially concerned about this reduced precision, as I seldom want my continued fractions to go out that far. And when I do, I have Mathematica to fall back on.

Finally, a small advantage of doing continued fractions in Python instead of Mathematica is that Python uses zero-based indexing for lists, which is consistent with the standard notation given above. Mathematica uses one-based indexing, which usually works out nicely when dealing with vectors and matrices but not in this case.

Update 16 Aug 2023 10:53 PM
Shortly after this post was published Thomas J. Fan got in touch with me on Mastodon and told me that the SymPy library has continued fraction functions in the number theory sublibrary. Of course! I felt silly for not looking there.

He also included this code snippet, which I ran in Jupyter:

In [1]: from itertools import islice

In [2]: from sympy.core import pi

In [3]: from sympy.ntheory.continued_fraction import continued_fraction_iterator

In [4]: list(islice(continued_fraction_iterator(pi), 15))
Out[4]: [3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1]

So the function is not only prewritten, it’s more accurate than mine. Of course, the final line to get a list of terms is kind of convoluted, but it could easily be wrapped in a more compact function. I may give it a go.

Thanks, Thomas!

Update 18 Aug 2023 9:49 AM
I played around with SymPy’s continued fraction functions yesterday and have decided to stick with my function. As discussed above, the SymPy functions are more accurate than mine, but to get that accuracy you have to be working with the SymPy definitions of numbers like $π$ and functions like sqrt. Since I’m usually working with the math package’s definitions, which are numeric rather than symbolic, I wouldn’t normally get the value out of the SymPy functions. Still, it’s good to know that they’re there.

It’s a thin Dover paperback, so it’s pretty cheap. ↩

And now it’s all this

I just said what I said and it was wrong
Or was taken wrong

Continued fractions in Python

Site search

Meta

Recent posts

Credits

And now it’s all this

I just said what I said and it was wrong Or was taken wrong

Continued fractions in Python

Site search

Meta

Recent posts

Credits

I just said what I said and it was wrong
Or was taken wrong