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) for _ in xrange(300)][5::5]
2 thoughts on “python: continued fractions”
Only quadratic irrationals have continued fraction expansions which are eventually periodic. Other irrational numbers of higher degree do not have period expansions.
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.