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


 
// Little Man Computer, for Rosetta Code.
// Displays terms of Calkin-Wilf sequence up to the given index.
// The chosen algorithm calculates the i-th term directly from i
//  (i.e. not using any previous terms).
input    INP        // get number of terms from user
         BRZ exit   // exit if 0
         STA max_i  // store maximum index
         LDA c1     // index := 1
next_i   STA i
// Write index followed by '->'
         OTX 3      // non-standard: minimum width 3, no new line
         LDA asc_hy
         OTC
         LDA asc_gt
         OTC
// Find greatest power of 2 not exceeding i,
//  and count the number of binary digits in i.
         LDA c1
         STA pwr2
loop2    STA nrDigits
         LDA i
         SUB pwr2
         SUB pwr2
         BRP double
         BRA part2    // jump out if next power of 2 would exceed i
double   LDA pwr2
         ADD pwr2
         STA pwr2
         LDA nrDigits
         ADD c1
         BRA loop2
// The nth term a/b is calculated from the binary digits of i.
// The leading 1 is not used.
part2    LDA c1
         STA a     // a := 1
         STA b     // b := 1
         LDA i
         SUB pwr2
         STA diff
// Pre-decrement count, since leading 1 is not used
dec_ct   LDA nrDigits  // count down the number of digits
         SUB c1
         BRZ output    // if all digits done, output the result
         STA nrDigits
// We now want to compare diff with pwr2/2.
// Since division is awkward in LMC, we compare 2*diff with pwr2.
         LDA diff      // diff := 2*diff
         ADD diff
         STA diff
         SUB pwr2      // is diff >= pwr2 ?
         BRP digit_1   // binary digit is 1 if yes, 0 if no
// If binary digit is 0 then set b := a + b
         LDA a
         ADD b
         STA b
         BRA dec_ct
// If binary digit is 1 then update diff and set a := a + b
digit_1  STA diff
         LDA a
         ADD b
         STA a
         BRA dec_ct
// Now have nth term a/b. Write it to the output.
output   LDA a         // write a
         OTX 1         // non-standard: minimum width 1; no new line
         LDA asc_sl    // write slash
         OTC
         LDA b         // write b
         OTX 11        // non-standard: minimum width 1; add new line
         LDA i         // have we done maximum i yet?
         SUB max_i
         BRZ exit      // if yes, exit
         LDA i         // if no, increment i and loop back
         ADD c1
         BRA next_i
exit     HLT
// Constants
c1       DAT 1
asc_hy   DAT 45
asc_gt   DAT 62
asc_sl   DAT 47
// Variables
i        DAT
max_i    DAT
pwr2     DAT
nrDigits DAT
diff     DAT
a        DAT
b        DAT
// end
 