Python 2.7 - Continued Fraction Expansion - Understanding the error

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?

aneroid · Accepted Answer · 2016-06-29T11:42:46.237

4

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

edited Jun 29 '16 at 11:42

answered Jun 29 '16 at 11:08

aneroid

12,983
3
36
66

I don't understand what you're staying, `N//1` is not equivalent to `N` i.e. 3.223//1 = 3. – ggordon Jun 29 '16 at 11:13
Again I don't see your point, if I use `int(N)` instead of `N//1` I receive the same error. With `N = 1/f` I am trying to get the reciprocal of f, not use floor division. If I print 1/f for the first 10 expansions I get: `4.08163265306 12.25 4.0 1.0 2.47777268231e+11 4.73320814676 1.36386918834 2.74824038982 1.33646888567 2.97204301075` – ggordon Jun 29 '16 at 11:22
Ah, ok I see your point now, had't yet seen your changes to your post – ggordon Jun 29 '16 at 11:24
I meant use `int(N)` instead of `N//1` only because it's clearer. Python will do float math when required. Although, appears less performant (will add that code). – aneroid Jun 29 '16 at 11:24
Testing for precision seems a bit messy, I thought there might be a neater way to do this, but thanks – ggordon Jun 29 '16 at 11:29
The joy of floats :-) And if this answers your question, feel free to Accept. – aneroid Jun 29 '16 at 11:52

Jerzyk · Answer 2 · 2016-06-29T11:44:28.770

you are using float for your operation, unfortunately some of the numbers can't be represented as a number in binary representation.

there are two options how to fix it, first - assume that your numbers are "close enough" (even new Python 3.5.2 introduces math.isclose), or you are using different implementation of floats e.g. Decimal you can get proper results.

N.b. this is why for all financial systems, nobody is ever using floats, only int/bigint or Decimals.

 In [21] > N = decimal.Decimal('3.245')

 In [22] > while True:
    print 'N: %s' % (N//1,)
    f = N - N//1
    print 'f: %s' % f
    if f == 0:
        break
    N = 1/f

N: 3
f: 0.245
N: 4
f: 0.081632653061224489795918367
N: 12
f: 0.25000000000000000000000005
N: 3
f: 0.999999999999999999999999200
N: 1
f: 8.00E-25
N: 1250000000000000000000000
f: 0

Good suggestion to use `math.isclose` but, as you pointed out, it's Python 3.x and the OP is using the future division, so he/she is on Py2k. — aneroid, Jun 29 '16 at 11:46

score -1 · Answer 3 · answered Jun 29 '16 at 11:32

Apparently the bug is due to the integer division of 4.0 by 1. 4.0 in floating point representation has a error involved (if you have idea about floating point representation). So actually int(4.0) < 4. This causes the N//1 to become 3.0 and the fraction f is something like 0.999999.... So this never ends. Print the N and f in every iteration and play around. You will realize.

Python 2.7 - Continued Fraction Expansion - Understanding the error

3 Answers3