Continued fraction convergents

Given a positive real number, if we truncate its continued fraction representation at a certain depth, we obtain a rational approximation to the real number. The sequence of successively better such approximations is its convergent sequence.

Continued fraction convergents is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Problem:

• Given a positive rational number ${\frac {m}{n}}$, specified by two positive integers $m,n$, output its entire sequence of convergents.
• Given a quadratic real number ${\frac {b{\sqrt {a}}+m}{n}}>0$, specified by integers $a,b,m,n$, where $a$ is not a perfect square, output the first $k$ convergents when given a positive number $k$.

The output format can be whatever is necessary to represent rational numbers, but it probably should be a 2-tuple of integers.

For example, given $a=2,b=1,m=0,n=1$, since

$${\sqrt {2}}=1+{\cfrac {1}{2+{\cfrac {1}{2+{\cfrac {1}{2+\ddots }}}}}}$$

the program should output $(1,1),(3,2),(7,5),(17/12),(41/29),\dots$.

A simple check is to do this for the golden ratio ${\frac {{\sqrt {5}}+1}{2}}$, that is, $a=5,b=1,m=1,n=2$, which should output $(1,1),(2,1),(3,2),(5,3),(8,5),\dots$.

Print the results for 415/93, 649/200, ${\sqrt {2}}$, ${\sqrt {5}}$, and the golden ratio.

References and related tasks

• Wikipedia: Continued fraction
• Continued fraction
• Continued fraction/Arithmetic
• Continued fraction/Arithmetic/Construct from rational number
• Calkin-Wilf sequence

Phix

Translation of: Wren

with javascript_semantics
requires("1.0.5") -- mpq_get_d() added
include mpfr.e

procedure cfcRat(integer m, n)
    sequence p = {mpq_init(0), mpq_init(1)},
             q = {mpq_init(1), mpq_init(0)},
             s = {{sprintf("= %d/%d =",{m,n}),m/n}}
    mpq r = mpq_init_set_si(m, n),
      rem = mpq_init_set(r)
    while true do
        mpq whole = mpq_init_set_si(trunc(mpq_get_d(rem)))
        mpq {pn, qn, sn} = mpq_inits(3)
        mpq_mul(pn,whole,p[-1])
        mpq_add(pn,pn,p[-2])
        mpq_mul(qn,whole,q[-1])
        mpq_add(qn,qn,q[-2])
        mpq_div(sn,pn,qn)
        p &= pn
        q &= qn
        s &= {{mpq_get_str(sn),mpq_get_d(sn)}}
        if mpq_cmp(r,sn)=0  then exit end if
        mpq_sub(rem,rem,whole)
        mpq_inv(rem,rem)
    end while
    printf(1,"%s\n",join(s,"\n",fmt:="%-11s = %f"))
end procedure

procedure cfcQuad(string d, atom v, integer a, b, m, n, k)
    sequence p = {0, 1},
             q = {1, 0},
             s = {{sprintf("= %s =",d),v}}
    atom r = (sqrt(a)*b + m) / n,
       rem = r
    for i=1 to k do
        integer whole = trunc(rem),
                pn = whole * p[-1] + p[-2],
                qn = whole * q[-1] + q[-2]
        mpq sn = mpq_init_set_si(pn, qn)
        p &= pn
        q &= qn
        s &= {{mpq_get_str(sn),mpq_get_d(sn)}}
        rem = 1/(rem-whole)
    end for
    printf(1,"%s\n",join(s,"\n",fmt:="%-11s = %f"))
end procedure

cfcRat(415,93)
cfcRat(649,200)
cfcQuad("sqrt(2)",sqrt(2),2, 1, 0, 1, 8)
cfcQuad("sqrt(5)",sqrt(5),5, 1, 0, 1, 8)
cfcQuad("phi",(sqrt(5)+1)/2,5, 1, 1, 2, 8)

Output:

= 415/93 =  = 4.462366
4           = 4.000000
9/2         = 4.500000
58/13       = 4.461538
415/93      = 4.462366
= 649/200 = = 3.245000
3           = 3.000000
13/4        = 3.250000
159/49      = 3.244898
649/200     = 3.245000
= sqrt(2) = = 1.414214
1           = 1.000000
3/2         = 1.500000
7/5         = 1.400000
17/12       = 1.416667
41/29       = 1.413793
99/70       = 1.414286
239/169     = 1.414201
577/408     = 1.414216
= sqrt(5) = = 2.236068
2           = 2.000000
9/4         = 2.250000
38/17       = 2.235294
161/72      = 2.236111
682/305     = 2.236066
2889/1292   = 2.236068
12238/5473  = 2.236068
51841/23184 = 2.236068
= phi =     = 1.618034
1           = 1.000000
2           = 2.000000
3/2         = 1.500000
5/3         = 1.666667
8/5         = 1.600000
13/8        = 1.625000
21/13       = 1.615385
34/21       = 1.619048

Wren

Library: Wren-rat

The following is loosely based on the Python code here. If a large number of terms were required for quadratic real numbers, then one might need to use 'arbitrary precision' arithmetic to minimize round-off errors when converting between floats and rationals.

import "./rat" for Rat

var cfcRat = Fn.new { |m, n|
    var p = [0, 1]
    var q = [1, 0]
    var s = []
    var r = Rat.new(m, n)
    var rem = r
    while (true) {
        var whole = rem.truncate
        var frac  = rem.fraction
        var pn = whole * p[-1] + p[-2]
        var qn = whole * q[-1] + q[-2]
        var sn = pn / qn
        p.add(pn)
        q.add(qn)
        s.add(sn)
        if (r == sn) break
        rem = frac.inverse
    }
    return s
}

var cfcQuad = Fn.new { |a, b, m, n, k|
    var p = [0, 1]
    var q = [1, 0]
    var s = []
    var rem = (a.sqrt * b + m) / n
    for (i in 1..k) {
        var whole = rem.truncate
        var frac  = rem.fraction
        var pn = whole * p[-1] + p[-2]
        var qn = whole * q[-1] + q[-2]
        var sn = Rat.new(pn, qn)
        p.add(pn)
        q.add(qn)
        s.add(sn)
        rem = 1 / frac
    }
    return s
}

System.print("The continued fraction convergents for the following (maximum 8 terms) are:")
System.print("415/93  = %(cfcRat.call(415, 93))")
System.print("649/200 = %(cfcRat.call(649, 200))")
System.print("√2      = %(cfcQuad.call(2, 1, 0, 1, 8))") 
System.print("√5      = %(cfcQuad.call(5, 1, 0, 1, 8))")
System.print("phi     = %(cfcQuad.call(5, 1, 1, 2, 8))")

Output:

The continued fraction convergents for the following (maximum 8 terms) are:
415/93  = [4/1, 9/2, 58/13, 415/93]
649/200 = [3/1, 13/4, 159/49, 649/200]
√2      = [1/1, 3/2, 7/5, 17/12, 41/29, 99/70, 239/169, 577/408]
√5      = [2/1, 9/4, 38/17, 161/72, 682/305, 2889/1292, 12238/5473, 51841/23184]
phi     = [1/1, 2/1, 3/2, 5/3, 8/5, 13/8, 21/13, 34/21]