python: continued fractions

When transferring pen-and-paper calculations into computer code, issues related to floating point precision tend to arise. For my preferred rapid prototyping language, python, this is the case as well. Of course, you can work with the decimal module, but this does not help you out in each and every case. For example, let's consider continued fractions.

Every (added, see comments: quadratic) irrational but not transcendental number can be expressed as a continued fraction which is partially periodical. The euclidean algorithm to calculate the corresponding expansion is simple but requires taking the inverse at every step and is, therefore, prone to numerical errors. In fact, even with arbitrary precision you do not work with exact digital numbers but with an infinite amount of numbers on a interval that map to the same digital representation. You can play around with the following code to see the problem: the last line should print only identical values (for sqrt(13)) but does not. Therefore, you have to test for periodicity with an tolerance limit although you work with (potentially) arbitrary precision numbers.

#!/usr/bin/env python
import decimal

def continued_fraction(n, precision=28):
    decimal.getcontext().prec = precision
    alpha = decimal.Decimal(n)

    while True:
        q, r = divmod(alpha, 1)
        yield (q, alpha)
        if r.is_zero():

        alpha = 1/(alpha-q)

def continued_fraction_float(n):
    alpha = n

    while True:
        q, r = divmod(alpha, 1)
        yield (q, alpha)
        if r == 0:

        alpha = 1/(alpha-q)

def detect_periodicity(iterable, limit=50, etol=0):
    states = []
    values = []

    is_periodic = False
    periodic_state = None
    while len(values) <= limit:
            value, state = next(iterable)
        for former in states:
            if abs(former-state) <= etol:
                is_periodic = True
                periodic_state = former
        if is_periodic:

    if is_periodic:
        split = states.index(periodic_state)
        return (is_periodic, values[:split], values[split:])
        return (is_periodic, values, [])

decimal.getcontext().prec = precision
d = decimal.Decimal(13)
d = d.sqrt()

g = continued_fraction(d)
is_periodic, non_periodic_elements, periodic_elements = detect_periodicity(g, etol=decimal.Decimal(1e-9))

g = continued_fraction(d, precision=100)
# should print only identical values - unfortunately, it does not
print [next(g)[1] for _ in xrange(300)][5::5] 

You may also like...

2 Responses

  1. R says:

    Only quadratic irrationals have continued fraction expansions which are eventually periodic. Other irrational numbers of higher degree do not have period expansions.

    • Guido says:

      Thank you for pointing this out. Apparently all quadratic irrationals have periodic continued fractions and periodic representations are only possible for quadratic irrationals. I have updated the post to make this more precise.

Leave a Reply

Your email address will not be published.