Continued fraction convergents

From Rosetta Code
Revision as of 11:26, 13 February 2024 by Petelomax (talk | contribs) (→‎{{header|Phix}}: made output more consistent with other entries)
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.

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.

Problem:

  • Given a positive rational number Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac mn} , specified by two positive integers Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle m, n} , output its entire sequence of convergents.
  • Given a quadratic real number Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{b\sqrt{a} + m}{n} > 0} , specified by integers Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a, b, m, n } , where Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} is not a perfect square, output the first Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle k} convergents when given a positive number Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a=2, b=1, m = 0, n = 1} , sinceFailed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sqrt{2} = 1 + \cfrac{1}{2 + \cfrac{1}{2 + \cfrac{1}{2 + \ddots}}}} the program should output Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (1, 1), (3,2), (7,5), (17/12), (41/29), \dots} .

A simple check is to do this for the golden ratio Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\sqrt{5}+1}{2}} , that is, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a=5, b=1, m = 1, n = 2 } , which should output Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (1,1), (2,1), (3,2), (5,3), (8,5), \dots} .

Print the results for 415/93, 649/200, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sqrt{2}} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sqrt{5}} , and the golden ratio.

References and related tasks


Julia 

function convergents(x::Real, maxcount::T) where T <: Integer
    components = T[]
    rationals = Rational{T}[]
    for _ in 1:maxcount
        fpart, ipart = modf(x)
        push!(components, T(ipart))
        fpart == 0 && break
        x = inv(fpart)
    end
    numa, denoma, numb, denomb = T(1), T(0), T(components[begin]), T(1)
    push!(rationals, numb // denomb)
    for comp in components[begin+1:end]
        numa, denoma, numb, denomb = numb, denomb, numa + comp * numb, denoma + comp * denomb
        push!(rationals, numb // denomb)
    end
    return rationals
end   

const tests = [("415/93", 415//93), ("649/200", 649//200), ("sqrt(2)", sqrt(2)),
               ("sqrt(5)", sqrt(5)), ("golden ratio", (sqrt(5) + 1) / 2)]

println("The continued fraction convergents for the following (maximum 8 terms) are:")
for (s, x) in tests
    println(lpad(s, 15), " = ", convergents(x, 8))
end
Output:
The continued fraction convergents for the following (maximum 8 terms) are:
         415/93 = Rational{Int64}[4, 9//2, 58//13, 415//93]
        649/200 = Rational{Int64}[3, 13//4, 159//49, 649//200]
        sqrt(2) = Rational{Int64}[1, 3//2, 7//5, 17//12, 41//29, 99//70, 239//169, 577//408]
        sqrt(5) = Rational{Int64}[2, 9//4, 38//17, 161//72, 682//305, 2889//1292, 12238//5473, 51841//23184]
   golden ratio = Rational{Int64}[1, 2, 3//2, 5//3, 8//5, 13//8, 21//13, 34//21]

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 = {},
             t = sprintf("%d/%d",{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)}
        if mpq_cmp(r,sn)=0  then exit end if
        mpq_sub(rem,rem,whole)
        mpq_inv(rem,rem)
    end while
    printf(1,"%14s = %s\n",{t,join(s)})
end procedure

procedure cfcQuad(string d, integer a, b, m, n, k)
    sequence p = {0, 1},
             q = {1, 0},
             s = {}
    atom rem = (sqrt(a)*b + m) / n
    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)}
        rem = 1/(rem-whole)
    end for
    printf(1,"%14s = %s\n",{d,join(s)})
end procedure

cfcRat(415,93)
cfcRat(649,200)
cfcQuad("sqrt(2)",2, 1, 0, 1, 8)
cfcQuad("sqrt(5)",5, 1, 0, 1, 8)
cfcQuad("golden ratio",5, 1, 1, 2, 8)
Output:
        415/93 = 4 9/2 58/13 415/93
       649/200 = 3 13/4 159/49 649/200
       sqrt(2) = 1 3/2 7/5 17/12 41/29 99/70 239/169 577/408
       sqrt(5) = 2 9/4 38/17 161/72 682/305 2889/1292 12238/5473 51841/23184
  golden ratio = 1 2 3/2 5/3 8/5 13/8 21/13 34/21

Raku

Translation of: Julia
# 20240210 Raku programming solution

sub convergents(Real $x is copy, Int $maxcount) {
   my @components = gather for ^$maxcount {
      my $fpart = $x - take $x.Int;
      $fpart == 0 ?? ( last ) !! ( $x = 1 / $fpart )
   }
   my ($numa, $denoma, $numb, $denomb) = 1, 0, @components[0], 1;
   return [ Rat.new($numb, $denomb) ].append: gather for @components[1..*] -> $comp {
      ( $numa, $denoma, $numb                , $denomb                  ) 
      = $numb, $denomb, $numa + $comp * $numb, $denoma + $comp * $denomb;
      take Rat.new($numb, $denomb);
   }
}

my @tests = [ "415/93", 415/93, "649/200", 649/200, "sqrt(2)", sqrt(2),
              "sqrt(5)", sqrt(5), "golden ratio", (sqrt(5) + 1) / 2     ];

say "The continued fraction convergents for the following (maximum 8 terms) are:";
for @tests -> $s, $x {
   say $s.fmt('%15s') ~ " = { convergents($x, 8).map: *.nude.join('/') } ";
}
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 
        sqrt(2) = 1/1 3/2 7/5 17/12 41/29 99/70 239/169 577/408 
        sqrt(5) = 2/1 9/4 38/17 161/72 682/305 2889/1292 12238/5473 51841/23184 
   golden ratio = 1/1 2/1 3/2 5/3 8/5 13/8 21/13 34/21

You may Attempt This Online!

Sidef 

func num2cfrac(n, r) {
    gather {
        r.times {
            n = 1/((n - take(n.floor.int)) || break)
        }
    }
}

func convergents(x, n) {
    var cfrac = num2cfrac(x, n)

    var(n1, n2) = (0, 1)
    var(d1, d2) = (1, 0)

    gather {
        for z in (cfrac) {
            (n1, n2) = (n2, n2*z + n1)
            (d1, d2) = (d2, d2*z + d1)
            take(n2/d2)
        }
    }
}

var tests = ["415/93", 415/93, "649/200", 649/200, "sqrt(2)", sqrt(2),
             "sqrt(5)", sqrt(5), "golden ratio", (sqrt(5) + 1) / 2   ]

var terms = 8
say "The continued fraction convergents for the following (maximum #{terms} terms) are:"
tests.each_slice(2, {|s,x|
    printf("%15s = %s\n", s, convergents(x, terms).map { .as_frac }.join(' '))
})
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
        sqrt(2) = 1/1 3/2 7/5 17/12 41/29 99/70 239/169 577/408
        sqrt(5) = 2/1 9/4 38/17 161/72 682/305 2889/1292 12238/5473 51841/23184
   golden ratio = 1/1 2/1 3/2 5/3 8/5 13/8 21/13 34/21

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]
Retrieved from "https://rosettacode.org/wiki/Continued_fraction_convergents?oldid=360717"