Python 2.7 - Continued Fraction Expansion - Understanding the error [duplicate]

Question

I've written this code to calculate the continued fraction expansion of a rational number N using the Euclidean algorithm:

from __future__ import division

def contFract(N):
    while True:
        yield N//1
        f = N - (N//1)
        if f == 0:
            break
        N = 1/f

If say N is 3.245 the function is never ending as apparently f never equals 0. The first 10 terms of the expansion are:

[3.0, 4.0, 12.0, 3.0, 1.0, 247777268231.0, 4.0, 1.0, 2.0, 1.0]

Which is clearly an error since the actual expansion is only:

[3;4,12,3,1] or [3;4,12,4]

What's causing the problem here? Is it some sort of rounding error?

Answer 1

The issue is you're testing f == 0 (integer 0) which is almost never true for floats. So the loop goes on forever.

Instead, check for some precision of equivalent to 0 (which can be wrong sometimes):

>>> from __future__ import division
>>>
>>> def contFract(N):
...     while True:
...         yield N//1
...         f = N - (N//1)
...         if f < 0.0001:  # or whatever precision you consider close enough to 0
...             break
...         N = 1/f
...
>>>
>>> list(contFract(3.245))
[3.0, 4.0, 12.0, 3.0, 1.0]
>>>

And in case f can be negative, do either -0.0001 < f < 0.0001 or abs(f) < 0.0001. Which is also considered inaccurate, see the linked article.

And wrt my comment to use int(N) instead of N//1 because it's clearer - it is slightly less efficient:

>>> import timeit
>>> timeit.timeit('N = 2.245; N//1', number=10000000)
1.5497028078715456
>>> timeit.timeit('N = 2.245; int(N)', number=10000000)
1.7633858824068103