Calkin-Wilf sequence: Difference between revisions

Content added Content deleted
(Added 3rd algorithm for calculating terms, and put all 3 algorithms into one program.)
Line 946: Line 946:
These programs were written in Free Pascal, using the Lazarus IDE and the Free Pascal compiler version 3.2.0. They are based on the Wikipedia article "Calkin-Wilf tree", rather than the algorithms in the task description.
These programs were written in Free Pascal, using the Lazarus IDE and the Free Pascal compiler version 3.2.0. They are based on the Wikipedia article "Calkin-Wilf tree", rather than the algorithms in the task description.
<lang pascal>
<lang pascal>
program CWTerms1;
program CWTerms;


{-------------------------------------------------------------------------------
{-------------------------------------------------------------------------------
Line 952: Line 952:
Calculates the Calkin-Wilf sequence up to the specified maximum index,
Calculates the Calkin-Wilf sequence up to the specified maximum index,
where the first term 1/1 has index 1.
where the first term 1/1 has index 1.
Command line format is
Command line format is: CWTerms <max_index>
CWTerms1 <max_index>
e.g. for the Rosetta Code example
CWTerms1 20


The program demonstrates 3 algorithms for calculating the sequence:
This program stores the terms in an array. Each term after the first
(1) Calculate term[2n] and term[2n + 1] from term[n]
is calculated from an earlier term in the array.
(2) Calculate term[n + 1] from term[n]
Cf. the alternative method in program CWTerms2.
(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.
-------------------------------------------------------------------------------}
-------------------------------------------------------------------------------}


Line 970: Line 971:
var
var
terms : array of TRational;
terms : array of TRational;
max_index, j, k : integer;
max_index, k : integer;

begin
// Read and validate maximum index
// Routine to calculate array of terms up the the maiximum index
procedure CalcTerms_algo_1();
max_index := SysUtils.StrToIntDef( paramStr(1), -1); // -1 if not an integer
var
if (max_index <= 0) then begin
WriteLn( 'Maximum index must be a positive integer');
j, k : integer;
exit;
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;
end;


// Method to get each term from the preceding term.
// Calculate the array of terms.
// A dynamic array such as "terms" has zero-based indices.
// a/b --> b/(a + b - 2(a mod b));
function CheckTerms_algo_2() : boolean;
// For convenience, we use indices [1..max_index] and ignore index [0].
var
SetLength( terms, max_index + 1);
j := 1; // index to earlier term, from which current term is calculated
index, a, b, temp : integer;
begin
k := 1; // index to current term
terms[1].Num := 1;
result := true;
terms[1].Den := 1;
index := 1;
a := 1;
while (k < max_index) do begin
inc(k);
b := 1;
while (index <= max_index) do begin
if (k and 1) = 0 then begin // or could write "if not Odd(k)"
terms[k].Num := terms[j].Num;
if (a <> terms[index].Num) or (b <> terms[index].Den) then
terms[k].Den := terms[j].Num + terms[j].Den;
result := false;
temp := a + b - 2*(a mod b);
end
else begin
a := b;
terms[k].Num := terms[j].Num + terms[j].Den;
b := temp;
inc( index)
terms[k].Den := terms[j].Den;
inc(j);
end;
end;
end;
end;


// Mathod to calcualte each term from its index, without using other terms.
// Display the terms
function CheckTerms_algo_3() : boolean;
for k := 1 to max_index do
var
with terms[k] do
index, a, b, pwr2, idiv2 : integer;
WriteLn( SysUtils.Format( '%8d: %d/%d', [k, Num, Den]));
begin
end.
result := true;
</lang>
for index := 1 to max_index do begin
{{out}}
Same as next program
<lang pascal>
program CWTerms2;


idiv2 := index div 2;
{-------------------------------------------------------------------------------
pwr2 := 1;
FreePascal command-line program.
while (pwr2 <= idiv2) do pwr2 := pwr2 shl 1;
Calculates the Calkin-Wilf sequence up to the specified maximum index,
where the first term 1/1 has index 1.
a := 1;
b := 1;
Command line format is
while (pwr2 > 1) do begin
CWTerms2 <max_index>
pwr2 := pwr2 shr 1;
e.g. for the Rosetta Code example
if (pwr2 and index) = 0 then
CWTerms2 20
inc( b, a)

else
This program calculates each term directly from its index,
inc( a, b);
without making use of previously calculated terms.
end;
Cf. the alternative method in program CWTerms1.
if (a <> terms[index].Num) or (b <> terms[index].Den) then
-------------------------------------------------------------------------------}
result := false;

end;
uses SysUtils;

// Subroutine to calculate the term at the passed-in index.
procedure CalculateTerm( index : integer; out num, den : integer);
var
a, b, pwr2, idiv2 : integer;
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;
end;
num := a;
den := b;
end;


// Main routine
var
max_index : integer;
k, num, den : integer;
begin
begin
// Read and validate maximum index
// Read and validate maximum index
Line 1,060: Line 1,050:
end;
end;


// Calculate and display the terms individually
// Calculate terms by algo 1, then check that algos 2 and 3 agree.
CalcTerms_algo_1();
for k := 1 to max_index do begin
if not CheckTerms_algo_2() then begin
CalculateTerm( k, {out} num, den);
WriteLn( SysUtils.Format( '%8d: %d/%d', [k, num, den]));
WriteLn( 'Algorithm 2 failed');
exit;
end;
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.
end.
</lang>
</lang>
Retrieved from "https://rosettacode.org/wiki/Calkin-Wilf_sequence"