Calkin-Wilf sequence: Difference between revisions
m (C++ - use a do/while loop) |
(this seems stable enough to be promoted from draft) |
||
| Line 1: | Line 1: | ||
{{ |
{{task}} |
||
The '''Calkin-Wilf sequence''' contains every nonnegative rational number exactly once. It can be calculated recursively as follows: |
The '''Calkin-Wilf sequence''' contains every nonnegative rational number exactly once. It can be calculated recursively as follows: |
||
Revision as of 08:46, 24 December 2020
You are encouraged to solve this task according to the task description, using any language you may know.
The Calkin-Wilf sequence contains every nonnegative rational number exactly once. It can be calculated recursively as follows:
a0 = 0
an+1 = 1/(2⌊an⌋+1-an) for n > 0
- Show on this page terms 0 through 20 of the Calkin-Wilf sequence. To avoid floating point error, you may want to use a rational number data type.
It is also possible, given a nonnegative rational number, to determine where it appears in the sequence without calculating the sequence. The procedure is to get the continued fraction representation of the rational and use it as the run-length encoding of the binary representation of the term number, beginning from the end of the continued fraction.
It only works if the number of terms in the continued fraction is odd- use either of the two equivalent representations to achieve this:
[a0; a1, a2, ..., an] = [a0; a1, a2 ,..., an-1, 1]
Thus, for example, the fraction 9/4 has odd continued fraction representation 2; 3, 1, giving a binary representation of 100011, which means 9/4 appears as the 35th term of the sequence.
- Find the position of the number 83116/51639 in the Calkin-Wilf sequence.
- See also
- Wikipedia entry: Calkin-Wilf tree
C++
<lang cpp>#include <iostream>
- include <vector>
- include <boost/rational.hpp>
using rational = boost::rational<unsigned long>;
unsigned long floor(const rational& r) {
return r.numerator()/r.denominator();
}
rational calkin_wilf_next(const rational& term) {
return 1UL/(2UL * floor(term) + 1UL - term);
}
std::vector<unsigned long> continued_fraction(const rational& r) {
unsigned long a = r.numerator();
unsigned long b = r.denominator();
std::vector<unsigned long> result;
do {
result.push_back(a/b);
unsigned long c = a;
a = b;
b = c % b;
} while (a != 1);
if (result.size() > 0 && result.size() % 2 == 0) {
--result.back();
result.push_back(1);
}
return result;
}
unsigned long term_number(const rational& r) {
unsigned long result = 0;
unsigned long d = 1;
unsigned long p = 0;
for (unsigned long n : continued_fraction(r)) {
for (unsigned long i = 0; i < n; ++i, ++p)
result |= (d << p);
d = !d;
}
return result;
}
int main() {
rational term = 0;
std::cout << "First 21 terms of the Calkin-Wilf sequence are:\n";
for (int i = 0; i <= 20; ++i) {
std::cout << std::setw(2) << i << ": " << term << '\n';
term = calkin_wilf_next(term);
}
rational r(83116, 51639);
std::cout << r << " is the " << term_number(r) << "th term of the sequence.\n";
}</lang>
- Output:
First 21 terms of the Calkin-Wilf sequence are: 0: 0/1 1: 1/1 2: 1/2 3: 2/1 4: 1/3 5: 3/2 6: 2/3 7: 3/1 8: 1/4 9: 4/3 10: 3/5 11: 5/2 12: 2/5 13: 5/3 14: 3/4 15: 4/1 16: 1/5 17: 5/4 18: 4/7 19: 7/3 20: 3/8 83116/51639 is the 123456789th term of the sequence.
Factor
<lang factor>USING: formatting io kernel lists lists.lazy math math.continued-fractions math.functions math.parser prettyprint sequences strings vectors ;
- next-cw ( x -- y ) [ floor dup + ] [ 1 swap - + recip ] bi ;
- calkin-wilf ( -- list ) 1 [ next-cw ] lfrom-by ;
- >continued-fraction ( x -- seq )
1vector [ dup last integer? ] [ dup next-approx ] until dup length even? [ unclip-last 1 - suffix! 1 suffix! ] when ;
- cw-index ( x -- n )
>continued-fraction <reversed> [ even? CHAR: 1 CHAR: 0 ? <string> ] map-index concat bin> ;
! Task "First 20 terms of the Calkin-Wilf sequence:" print 20 calkin-wilf ltake [ pprint bl ] leach nl nl
83116/51639 cw-index "83116/51639 is at index %d.\n" printf</lang>
- Output:
First 20 terms of the Calkin-Wilf sequence: 1 1/2 2 1/3 1+1/2 2/3 3 1/4 1+1/3 3/5 2+1/2 2/5 1+2/3 3/4 4 1/5 1+1/4 4/7 2+1/3 3/8 83116/51639 is at index 123456789.
FreeBASIC
Uses the code from Greatest common divisor#FreeBASIC as an include.
<lang freebasic>#include "gcd.bas"
type rational
num as integer den as integer
end type
dim shared as rational ONE, TWO ONE.num = 1 : ONE.den = 1 TWO.num = 2 : TWO.den = 1
function simplify( byval a as rational ) as rational
dim as uinteger g = gcd( a.num, a.den )
a.num /= g : a.den /= g
if a.den < 0 then
a.den = -a.den
a.num = -a.num
end if
return a
end function
operator + ( a as rational, b as rational ) as rational
dim as rational ret ret.num = a.num * b.den + b.num*a.den ret.den = a.den * b.den return simplify(ret)
end operator
operator - ( a as rational, b as rational ) as rational
dim as rational ret ret.num = a.num * b.den - b.num*a.den ret.den = a.den * b.den return simplify(ret)
end operator
operator * ( a as rational, b as rational ) as rational
dim as rational ret ret.num = a.num * b.num ret.den = a.den * b.den return simplify(ret)
end operator
operator / ( a as rational, b as rational ) as rational
dim as rational ret ret.num = a.num * b.den ret.den = a.den * b.num return simplify(ret)
end operator
function floor( a as rational ) as rational
dim as rational ret ret.den = 1 ret.num = a.num \ a.den return ret
end function
function cw_nextterm( q as rational ) as rational
dim as rational ret = (TWO*floor(q)) ret = ret + ONE : ret = ret - q return ONE / ret
end function
function frac_to_int( byval a as rational ) as uinteger
redim as uinteger cfrac(-1)
dim as integer lt = -1, ones = 1, ret = 0
do
lt += 1
redim preserve as uinteger cfrac(0 to lt)
cfrac(lt) = floor(a).num
a = a - floor(a) : a = ONE / a
loop until a.num = 0 or a.den = 0
if lt mod 2 = 1 and cfrac(lt) = 1 then
lt -= 1
cfrac(lt)+=1
redim preserve as uinteger cfrac(0 to lt)
end if
if lt mod 2 = 1 and cfrac(lt) > 1 then
cfrac(lt) -= 1
lt += 1
redim preserve as uinteger cfrac(0 to lt)
cfrac(lt) = 1
end if
for i as integer = lt to 0 step -1
for j as integer = 1 to cfrac(i)
ret *= 2
if ones = 1 then ret += 1
next j
ones = 1 - ones
next i
return ret
end function
function disp_rational( a as rational ) as string
if a.den = 1 or a.num= 0 then return str(a.num) return str(a.num)+"/"+str(a.den)
end function
dim as rational q q.num = 0 q.den = 1 for i as integer = 0 to 20
print i, disp_rational(q) q = cw_nextterm(q)
next i
q.num = 83116 q.den = 51639 print disp_rational(q)+" is the "+str(frac_to_int(q))+"th term."</lang>
- Output:
0 0 1 1 2 1/2 3 2 4 1/3 5 3/2 6 2/3 7 3 8 1/4 9 4/3 10 3/5 11 5/2 12 2/5 13 5/3 14 3/4 15 4 16 1/5 17 5/4 18 4/7 19 7/3 20 3/8 83116/51639 is the 123456789th term.
Go
Go just has arbitrary precision rational numbers which we use here whilst assuming the numbers needed for this task can be represented exactly by the 64 bit built-in types. <lang go>package main
import (
"fmt" "math" "math/big" "strconv" "strings"
)
func calkinWilf(n int) []*big.Rat {
cw := make([]*big.Rat, n+1)
cw[0] = new(big.Rat)
one := big.NewRat(1, 1)
two := big.NewRat(2, 1)
for i := 1; i <= n; i++ {
t := new(big.Rat).Set(cw[i-1])
f, _ := t.Float64()
f = math.Floor(f)
t.SetFloat64(f)
t.Mul(t, two)
t.Sub(t, cw[i-1])
t.Add(t, one)
t.Inv(t)
cw[i] = new(big.Rat).Set(t)
}
return cw
}
func toContinued(r *big.Rat) []int {
a := r.Num().Int64()
b := r.Denom().Int64()
var res []int
for {
res = append(res, int(a/b))
t := a % b
a, b = b, t
if a == 1 {
break
}
}
return res
}
func getTermNumber(cf []int) int {
b := ""
d := "1"
for _, n := range cf {
b = strings.Repeat(d, n) + b
if d == "1" {
d = "0"
} else {
d = "1"
}
}
i, _ := strconv.ParseInt(b, 2, 64)
return int(i)
}
func commatize(n int) string {
s := fmt.Sprintf("%d", n)
if n < 0 {
s = s[1:]
}
le := len(s)
for i := le - 3; i >= 1; i -= 3 {
s = s[0:i] + "," + s[i:]
}
if n >= 0 {
return s
}
return "-" + s
}
func main() {
cw := calkinWilf(20)
fmt.Println("The first 21 terms of the Calkin-Wilf sequence are:")
for i := 0; i <= 20; i++ {
fmt.Printf("%2d: %s\n", i, cw[i].RatString())
}
fmt.Println()
r := big.NewRat(83116, 51639)
cf := toContinued(r)
tn := getTermNumber(cf)
fmt.Printf("%s is the %sth term of the sequence.\n", r.RatString(), commatize(tn))
}</lang>
- Output:
The first 21 terms of the Calkin-Wilf sequence are: 0: 0 1: 1 2: 1/2 3: 2 4: 1/3 5: 3/2 6: 2/3 7: 3 8: 1/4 9: 4/3 10: 3/5 11: 5/2 12: 2/5 13: 5/3 14: 3/4 15: 4 16: 1/5 17: 5/4 18: 4/7 19: 7/3 20: 3/8 83116/51639 is the 123,456,789th term of the sequence.
Julia
<lang julia>function calkin_wilf(n)
cw = zeros(Rational, n + 1)
for i in 2:n + 1
t = Int(floor(cw[i - 1])) * 2 - cw[i - 1] + 1
cw[i] = 1 // t
end
return cw
end
function continued(r::Rational)
a, b = r.num, r.den
res = []
while true
push!(res, Int(floor(a / b)))
a, b = b, a % b
a == 1 && break
end
return res
end
function term_number(cf)
b, d = "", "1"
for n in cf
b = d^n * b
d = (d == "1") ? "0" : "1"
end
return parse(Int, b, base=2)
end
const cw = calkin_wilf(20) println("The first 21 terms of the Calkin-Wilf sequence are: $cw")
const r = 83116 // 51639 const cf = continued(r) const tn = term_number(cf) println("$r is the $tn-th term of the sequence.")
</lang>
- Output:
The first 21 terms of the Calkin-Wilf sequence are: Rational[0//1, 1//1, 1//2, 2//1, 1//3, 3//2, 2//3, 3//1, 1//4, 4//3, 3//5, 5//2, 2//5, 5//3, 3//4, 4//1, 1//5, 5//4, 4//7, 7//3, 3//8] 83116//51639 is the 123456789-th term of the sequence.
Phix
<lang Phix>function calkin_wilf(integer len)
sequence cw = repeat(0,len)
integer n=0, d=1
for i=1 to len do
{n,d} = {d,(floor(n/d)*2+1)*d-n}
cw[i] = {n,d}
end for
return cw
end function
function to_continued_fraction(sequence r)
integer {a,b} = r
sequence res = {}
while true do
res &= floor(a/b)
{a, b} = {b, remainder(a,b)}
if a=1 then exit end if
end while
return res
end function
function get_term_number(sequence cf)
sequence b = {}
integer d = 1
for i=1 to length(cf) do
b &= repeat(d,cf[i])
d = 1-d
end for
return bits_to_int(b)
end function
sequence cw = calkin_wilf(20) printf(1,"The first 21 terms of the Calkin-Wilf sequence are:\n 0: 0\n") for i=1 to 20 do
integer {n,d} = cw[i]
string s = iff(d=1?sprintf("%d",n):sprintf("%d/%d",{n,d}))
printf(1,"%2d: %s\n", {i, s})
end for printf(1,"\n") sequence r = {83116,51639} sequence cf = to_continued_fraction(r) integer tn = get_term_number(cf) printf(1,"%d/%d is the %,d%s term of the sequence.\n", r&{tn,ord(tn)})</lang>
- Output:
The first 21 terms of the Calkin-Wilf sequence are: 0: 0 1: 1 2: 1/2 3: 2 4: 1/3 5: 3/2 6: 2/3 7: 3 8: 1/4 9: 4/3 10: 3/5 11: 5/2 12: 2/5 13: 5/3 14: 3/4 15: 4 16: 1/5 17: 5/4 18: 4/7 19: 7/3 20: 3/8 83116/51639 is the 123,456,789th term of the sequence.
Raku
Technically, the Calkin-Wilf sequence should begin with 1, but start with 0 as that is what the task specifies.
Conveniently, having the bogus first term shifts the indices up by one, making the ordinal position and index match.
Only show the first twenty terms that are actually in the sequence.
<lang perl6>my @calkin-wilf = 0, 1, {1 / (.Int × 2 + 1 - $_)} … *;
- Rational to Calkin-Wilf index
sub r2cw (Rat $rat) { :2( join , flat (flat (1,0) xx *) Zxx reverse r2cf $rat ) }
- The task
say "First twenty terms of the Calkin-Wilf sequence: ",
@calkin-wilf[1..20]».&prat.join: ', ';
say "\n99991st through 100000th: ",
(my @tests = @calkin-wilf[99_991 .. 100_000])».&prat.join: ', ';
say "\nCheck reversibility: ", @tests».Rat».&r2cw.join: ', ';
say "\n83116/51639 is at index: ", r2cw 83116/51639;
- Helper subs
sub r2cf (Rat $rat is copy) { # Rational to continued fraction
gather loop {
$rat -= take $rat.floor; last if !$rat; $rat = 1 / $rat;
}
}
sub prat ($num) { # pretty Rat
return $num unless $num ~~ Rat|FatRat; return $num.numerator if $num.denominator == 1; $num.nude.join: '/';
}</lang>
- Output:
First twenty terms of the Calkin-Wilf sequence: 1, 1/2, 2, 1/3, 3/2, 2/3, 3, 1/4, 4/3, 3/5, 5/2, 2/5, 5/3, 3/4, 4, 1/5, 5/4, 4/7, 7/3, 3/8 99991st through 100000th: 1085/303, 303/1036, 1036/733, 733/1163, 1163/430, 430/987, 987/557, 557/684, 684/127, 127/713 Check reversibility: 99991, 99992, 99993, 99994, 99995, 99996, 99997, 99998, 99999, 100000 83116/51639 is at index: 123456789
Rust
<lang rust>// [dependencies] // num = "0.3"
use num::rational::Rational;
fn calkin_wilf_next(term: &Rational) -> Rational {
Rational::from_integer(1) / (Rational::from_integer(2) * term.floor() + 1 - term)
}
fn continued_fraction(r: &Rational) -> Vec<isize> {
let mut a = *r.numer();
let mut b = *r.denom();
let mut result = Vec::new();
loop {
let (q, r) = num::integer::div_rem(a, b);
result.push(q);
a = b;
b = r;
if a == 1 {
break;
}
}
let len = result.len();
if len != 0 && len % 2 == 0 {
result[len - 1] -= 1;
result.push(1);
}
result
}
fn term_number(r: &Rational) -> usize {
let mut result: usize = 0;
let mut d: usize = 1;
let mut p: usize = 0;
for n in continued_fraction(r) {
for _ in 0..n {
result |= d << p;
p += 1;
}
d ^= 1;
}
result
}
fn main() {
println!("First 21 terms of the Calkin-Wilf sequence are:");
let mut term = Rational::from_integer(0);
for i in 0..=20 {
println!("{:2}: {}", i, term);
term = calkin_wilf_next(&term);
}
let r = Rational::new(83116, 51639);
println!("{} is the {}th term of the sequence.", r, term_number(&r));
}</lang>
- Output:
First 21 terms of the Calkin-Wilf sequence are: 0: 0 1: 1 2: 1/2 3: 2 4: 1/3 5: 3/2 6: 2/3 7: 3 8: 1/4 9: 4/3 10: 3/5 11: 5/2 12: 2/5 13: 5/3 14: 3/4 15: 4 16: 1/5 17: 5/4 18: 4/7 19: 7/3 20: 3/8 83116/51639 is the 123456789th term of the sequence.
Wren
<lang ecmascript>import "/rat" for Rat import "/fmt" for Fmt, Conv
var calkinWilf = Fn.new { |n|
var cw = List.filled(n + 1, null)
cw[0] = Rat.zero
for (i in 1..n) {
var t = cw[i-1].floor * 2 - cw[i-1] + 1
cw[i] = Rat.one / t
}
return cw
}
var toContinued = Fn.new { |r|
var a = r.num
var b = r.den
var res = []
while (true) {
res.add((a/b).floor)
var t = a % b
a = b
b = t
if (a == 1) break
}
return res
}
var getTermNumber = Fn.new { |cf|
var b = ""
var d = "1"
for (n in cf) {
b = (d * n) + b
d = (d == "1") ? "0" : "1"
}
return Conv.atoi(b, 2)
}
var cw = calkinWilf.call(20) System.print("The first 21 terms of the Calkin-Wilf sequence are:") Rat.showAsInt = true for (i in 0..20) Fmt.print("$2d: $s", i, cw[i]) System.print() var r = Rat.new(83116, 51639) var cf = toContinued.call(r) var tn = getTermNumber.call(cf) Fmt.print("$s is the $,r term of the sequence.", r, tn)</lang>
- Output:
The first 21 terms of the Calkin-Wilf sequence are: 0: 0 1: 1 2: 1/2 3: 2 4: 1/3 5: 3/2 6: 2/3 7: 3 8: 1/4 9: 4/3 10: 3/5 11: 5/2 12: 2/5 13: 5/3 14: 3/4 15: 4 16: 1/5 17: 5/4 18: 4/7 19: 7/3 20: 3/8 83116/51639 is the 123,456,789th term of the sequence.