The Calkin-Wilf sequence contains every nonnegative rational number exactly once.

It can be calculated recursively as follows:

       a1   =  1 
       an+1  =  1/(2⌊an⌋+1-an) for n > 1 



Task part 1

 Show on this page terms 1 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 non-negative 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] 



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.



Task part 2

 Find the position of the number   83116/51639   in the Calkin-Wilf sequence.


See also

 Wikipedia entry: Calkin-Wilf tree
 Continued fraction
 Continued fraction/Arithmetic/Construct from rational number


 
program CWTerms;
 
{-------------------------------------------------------------------------------
FreePascal command-line program.
Calculates the Calkin-Wilf sequence up to the specified maximum index,
  where the first term 1/1 has index 1.
Command line format is: CWTerms <max_index>
 
The program demonstrates 3 algorithms for calculating the sequence:
(1) Calculate term[2n] and term[2n + 1] from term[n]
(2) Calculate term[n + 1] from term[n]
(3) Calculate term[n] directly from n, without using other terms
Algorithm 1 is called first, and stores the terms in an array.
Then the program calls Algorithms 2 and 3, and checks that they agree
  with Algorithm 1.
-------------------------------------------------------------------------------}
 
uses SysUtils;
 
type TRational = record
  Num, Den : integer;
end;
 
var
  terms : array of TRational;
  max_index, k : integer;
 
  // Routine to calculate array of terms up the the maiximum index
  procedure CalcTerms_algo_1();
  var
    j, k : integer;
  begin
    SetLength( terms, max_index + 1);
    j := 1; // index to earlier term, from which current term is calculated
    k := 1; // index to current term
    terms[1].Num := 1;
    terms[1].Den := 1;
    while (k < max_index) do begin
      inc(k);
      if (k and 1) = 0 then begin // or could write "if not Odd(k)"
        terms[k].Num := terms[j].Num;
        terms[k].Den := terms[j].Num + terms[j].Den;
      end
      else begin
        terms[k].Num := terms[j].Num + terms[j].Den;
        terms[k].Den := terms[j].Den;
        inc(j);
      end;
    end;
  end;
 
  // Method to get each term from the preceding term.
  // a/b --> b/(a + b - 2(a mod b));
  function CheckTerms_algo_2() : boolean;
  var
    index, a, b, temp : integer;
  begin
    result := true;
    index := 1;
    a := 1;
    b := 1;
    while (index <= max_index) do begin
      if (a <> terms[index].Num) or (b <> terms[index].Den) then
        result := false;
      temp := a + b - 2*(a mod b);
      a := b;
      b := temp;
      inc( index)
    end;
  end;
 
  // Mathod to calcualte each term from its index, without using other terms.
  function CheckTerms_algo_3() : boolean;
  var
    index, a, b, pwr2, idiv2 : integer;
  begin
    result := true;
    for index := 1 to max_index do begin
 
      idiv2 := index div 2;
      pwr2 := 1;
      while (pwr2 <= idiv2) do pwr2 := pwr2 shl 1;
      a := 1;
      b := 1;
      while (pwr2 > 1) do begin
        pwr2 := pwr2 shr 1;
        if (pwr2 and index) = 0 then
          inc( b, a)
        else
          inc( a, b);
      end;
      if (a <> terms[index].Num) or (b <> terms[index].Den) then
        result := false;
    end;
  end;
 
begin
  // Read and validate maximum index
  max_index := SysUtils.StrToIntDef( paramStr(1), -1); // -1 if not an integer
  if (max_index <= 0) then begin
    WriteLn( 'Maximum index must be a positive integer');
    exit;
  end;
 
  // Calculate terms by algo 1, then check that algos 2 and 3 agree.
  CalcTerms_algo_1();
  if not CheckTerms_algo_2() then begin
    WriteLn( 'Algorithm 2 failed');
    exit;
  end;
  if not CheckTerms_algo_3() then begin
    WriteLn( 'Algorithm 3 failed');
    exit;
  end;
 
  // Display the terms
  for k := 1 to max_index do
    with terms[k] do
      WriteLn( SysUtils.Format( '%8d: %d/%d', [k, Num, Den]));
end.
 