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?
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