Source code in the pascal programming language

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.


program CWIndex;

{-------------------------------------------------------------------------------
FreePascal command-line program.
Calculates index of a rational number in the Calkin-Wilf sequence,
  where the first term 1/1 has index 1.
Command line format is
  CWIndex <numerator> <denominator>
e.g. for the Rosetta Code example
  CWIndex 83116 51639
-------------------------------------------------------------------------------}

uses SysUtils;

var
  num, den : integer;
  a, b : integer;
  pwr2, index : qword; // 64-bit unsiged
begin
  // Read and validate input.
  num := SysUtils.StrToIntDef( paramStr(1), -1); // return -1 if not an integer
  den := SysUtils.StrToIntDef( paramStr(2), -1);
  if (num <= 0) or (den <= 0) then begin
    WriteLn( 'Numerator and denominator must be positive integers');
    exit;
  end;

  // Input OK, calculate and display index of num/den
  // The index may overflow 64 bits, so turn on overflow detection
{$Q+}
  a := num;
  b := den;
  pwr2 := 1;
  index := 0;
  try
    while (a <> b) do begin
      if (a < b) then
        dec( b, a)
      else begin
        dec( a, b);
        inc( index, pwr2);
      end;
      pwr2 := 2*pwr2;
    end;
    inc( index, pwr2);
    WriteLn( SysUtils.Format( 'Index of %d/%d is %u', [num, den, index]));
  except
    WriteLn( 'Index is too large for 64 bits');
  end;
end.