convergents_to_cube_root_of_2.sf

#!/usr/bin/ruby

# Compute the convergents to cube root of 2.

# Formula due to Enrico Bombieri and Alfred J. van der Poorten:
#   https://web.williams.edu/Mathematics/sjmiller/public_html/book/papers/vdp/BombieriPoorten_CFofAlgNumbs.pdf

# Algorithm adapted after Robert Israel, Oct 08 2017.

# See also:
#   https://oeis.org/A002352 -- Numerators of convergents to cube root of 2.
#   https://oeis.org/A002351 -- Denominators of convergents to cube root of 2.

var (a, c, q)

c = -> (n) is cached {

    return 1 if (n == 0)

    floor((-1)**n * 3*a(n)**2 / (q(n)*(a(n)**3 - 2*q(n)**3)) - q(n-1)/q(n))
}

a = -> (n) is cached {

    return 1 if (n == 0)
    return 1 if (n == 1)

    c(n-1)*a(n-1) + a(n-2)
}

q = -> (n) is cached {

    return 0 if (n == 0)
    return 1 if (n == 1)

    c(n-1)*q(n-1) + q(n-2)
}

say ("Numerators  : ", 15.of { a(_) })
say ("Denominators: ", 15.of { q(_) })

say ''

for n in (1..20) {
    printf("2^(1/3) =~ %22s =~ %s\n", as_rat(a(n)/q(n)) a(n)/q(n))
}

__END__
Numerators  : [1, 1, 4, 5, 29, 34, 63, 286, 349, 635, 5429, 6064, 90325, 96389, 1054215]
Denominators: [0, 1, 3, 4, 23, 27, 50, 227, 277, 504, 4309, 4813, 71691, 76504, 836731]

2^(1/3) =~                      1 =~ 1
2^(1/3) =~                    4/3 =~ 1.33333333333333333333333333333333333333333333333
2^(1/3) =~                    5/4 =~ 1.25
2^(1/3) =~                  29/23 =~ 1.26086956521739130434782608695652173913043478261
2^(1/3) =~                  34/27 =~ 1.25925925925925925925925925925925925925925925926
2^(1/3) =~                  63/50 =~ 1.26
2^(1/3) =~                286/227 =~ 1.25991189427312775330396475770925110132158590308
2^(1/3) =~                349/277 =~ 1.25992779783393501805054151624548736462093862816
2^(1/3) =~                635/504 =~ 1.25992063492063492063492063492063492063492063492
2^(1/3) =~              5429/4309 =~ 1.25992109538175910884195869111162682757020190299
2^(1/3) =~              6064/4813 =~ 1.25992104716393102015375025971327654269686266362
2^(1/3) =~            90325/71691 =~ 1.2599210500620719476642814300260841667712823088
2^(1/3) =~            96389/76504 =~ 1.2599210498797448499424866673637979713479033776
2^(1/3) =~         1054215/836731 =~ 1.25992104989536661125260089562834411537280201164
2^(1/3) =~        2204819/1749966 =~ 1.25992104989468366813983814542682543546560333172
2^(1/3) =~        3259034/2586697 =~ 1.25992104989490458294883397630259748242643030861
2^(1/3) =~      15240955/12096754 =~ 1.25992104989487262450736784429938808377850785426
2^(1/3) =~    186150494/147747745 =~ 1.25992104989487318402050738574724101542125059168
2^(1/3) =~    387541943/307592244 =~ 1.25992104989487316201639986735166183188936324415
2^(1/3) =~  1348776323/1070524477 =~ 1.2599210498948731650532825696427303642119338482