implement continued fractions

[GoogleCodeExporter]

It'd be nice to have a way to convert any mpf number to a continued
fraction representation.

Original issue reported on code.google.com by ondrej.c...@gmail.com on 27 Aug 2008 at 11:47

[GoogleCodeExporter]

This should go via mpf -> mpq -> continued fraction.

The algorithm for converting an ordinary fraction to a continued fraction is 
very
elegant (like Euclid's algorithm), especially when rewritten to only use 
integers.
This is the shortest version I've come up with:

def contfrac(p, q):
    while q:
        n = p // q
        yield n
        q, p = p - q*n, q

>>> list(contfrac(225,100))
[2, 4]
>>> list(contfrac(-42,10))
[-5, 1, 4]
>>> list(contfrac(31415926536, 10000000000))
[3L, 7L, 15L, 1L, 292L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 45L, 1L, 1L, 8L]

Original comment by fredrik....@gmail.com on 27 Aug 2008 at 1:02

[GoogleCodeExporter]

That's really nice and easy.

Original comment by ondrej.c...@gmail.com on 27 Aug 2008 at 1:13

[GoogleCodeExporter]

def convergents(cf):
    r, s, p, q = 0, 1, 1, 0
    for c in cf:
        r, s, p, q = p, q, c*p+r, c*q+s
        yield p, q

>>> list(convergents(contfrac(31415926,10000000)))[:4]
[(3, 1), (22, 7), (333, 106), (355, 113)]

Unfortunately, this gives a "truncating" sequence of convergents; an additional 
trick
is required to obtain the true best approximation.

Original comment by fredrik....@gmail.com on 27 Aug 2008 at 1:34

[GoogleCodeExporter]

By the way, the truncating version works with symbolic expressions too:

>>> from sympy import *
>>> var('a b c d')
(a, b, c, d)
>>> pprint([p/q for (p, q) in convergents([a, b, c, d])])
    1 + a*b  a + c*(1 + a*b)  1 + a*b + d*(a + c*(1 + a*b))
[a, -------, ---------------, -----------------------------]
       b         1 + b*c             b + d*(1 + b*c)

Original comment by fredrik....@gmail.com on 27 Aug 2008 at 1:37

[GoogleCodeExporter]

Interesting. Go ahead with this, it will be very useful to have this.

Original comment by ondrej.c...@gmail.com on 27 Aug 2008 at 1:43

[GoogleCodeExporter]

Actually, forming best approximations from the "truncated" convergents *does* 
work,
if the first term is skipped. At least, I think so (but right now I am 
confused).
Here is a version that works with SymPy Rationals.

from sympy import *
from math import frexp

def contfrac(p, q):
    while q:
        n = p // q
        yield n
        q, p = p - q*n, q

def convergents(cf):
    p, q, r, s = 1, 0, 0, 1
    for c in cf:
        p, q, r, s = c*p+r, c*q+s, p, q
        yield p, q

def simplify(x, d=100):
    x, y = frexp(x)
    x = Integer(int(x*2**53)) * Integer(2)**(y-53)
    p, q = 0, 0
    for r, s in convergents(contfrac(x.p, x.q)):
        if s > d and q:
            break
        p, q = r, s
    return Rational(p, q)

simplify(0.1)
simplify(0.501)
simplify(0.499)
simplify(0.666)
simplify(7.2323)
simplify(3.1415926535897932)
simplify(3.1415926535897932, 1)
simplify(3.1415926535897932, 1000)

(d is the maximum acceptable denominator)

Original comment by fredrik....@gmail.com on 27 Aug 2008 at 4:55

[GoogleCodeExporter]

For the record, here is the output of it:

1/10
1/2
1/2
2/3
716/99
22/7
3
355/113

In my opinion, with the precision of floats in Python, I think one can get 
pretty far
with this.

Original comment by ondrej.c...@gmail.com on 27 Aug 2008 at 5:18

[GoogleCodeExporter]

For example if you try this:

a = float(53)/2420933
print a

print simplify(a, 1)
print simplify(a, 10)
print simplify(a, 100)
print simplify(a, 1000)
print simplify(a, 10000)
print simplify(a, 100000)
print simplify(a, 1000000)
print simplify(a, 10000000)
print simplify(a, 100000000)
print simplify(a, 10**14)
print simplify(a, 10**15)

it will output:

2.18923861173e-05
0
0
0
0
0
1/45678
1/45678
53/2420933
53/2420933
53/2420933
18591760240/849234073454791

BTW, 53/(2420933+1) = 1/45678. So we need some way how to choose the best 
answer. 

E.g. we can run some tests, e.g. try most of the all possible fractions, let's 
say
all fractions a/b with both a, b from -10**n to 10**n, where n = 5 or something 
that
the user sets. And we'll tune the algorithm to return the correct answer for all
these fractions. If you know what I mean.

Original comment by ondrej.c...@gmail.com on 27 Aug 2008 at 5:33

[GoogleCodeExporter]

Any progress on this? Let's get it in, seems to be working in the easy cases 
already.

Original comment by ondrej.c...@gmail.com on 4 Sep 2008 at 12:46

[GoogleCodeExporter]

[deleted comment]

[GoogleCodeExporter]

[deleted comment]

[GoogleCodeExporter]

[deleted comment]

[GoogleCodeExporter]

In the minds of most people continued fractions are tied up with integers and 
rationals. Not so for Euclid whose algorithm works perfectly well with floating
point numbers. Actually he might might have been thinking of ratio's of line 
segments, whit just fitting the smaller an integer number of times onto the 
larger. This amounts to repeatedly splitting x into floor(x) and a remainder.
If the remainder is zero the process ends.
x1 = floor(x) 
x = 1/(x-x1)
x2 = floor(x)
x= 1/(x-x2)
Bottom line x = x1 + 1/(x2 + 1/(x3 _ ...)

A brilliant reference about cd is gopher.txt (google for it, it is everywhere).

Original comment by alberhem...@gmail.com on 2 Jun 2014 at 10:07

[GoogleCodeExporter]

While I'm at it, I might as well show how its done without rationals:
-----------------
from sympy import *
from math import frexp

def contfrac( x):
    while True:
        n = floor(x)
        yield int(n)
        if n==x: return
        x=1/(x-n)

def convergents(cf):
    p, q, r, s = 1, 0, 0, 1
    for c in cf:
        p, q, r, s = c*p+r, c*q+s, p, q
        yield p, q

def simplify(x, d=100):
    p, q = 0, 0
    for r, s in convergents(contfrac(x)):
        if s > d and q:
            break
        p, q = r, s
    return Rational(p, q)

print    simplify(0.1)
print        simplify(0.501)
print        simplify(0.499)
print        simplify(0.666)
print        simplify(7.2323)
print        simplify(3.1415926535897932)
print        simplify(3.1415926535897932, 1)
print        simplify(3.1415926535897932, 3000)

# (d is the maximum acceptable denominator)
--------------------------------------
This gives the exact same output as in message #6

Original comment by alberhem...@gmail.com on 2 Jun 2014 at 8:03