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():
            return

        alpha = 1/(alpha-q)

def continued_fraction_float(n):
    alpha = n

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

        alpha = 1/(alpha-q)

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

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

    if is_periodic:
        split = states.index(periodic_state)
        return (is_periodic, values[:split], values[split:])
    else:
        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]