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


 
[For Rosetta Code. EDSAC program, Initial Orders 2.
 Prints the first 20 terms of the Calkin-Wilf sequence.
 Uses term{n} to calculate term{n + 1}.]
 
[Print subroutine for non-negative 17-bit integers.
 Parameters: 0F = integer to be printed (not preserved)
             1F = character for leading zero (preserved)
 Workspace: 4F, 5F. Even address; 40 locations]
            T   56 K [define load address]
GKA3FT34@A1FT35@S37@T36@AFT5FT4FH38#@V4DH30@
R32FR16FYFE23@O35@A2FT36@T5FV4DYFL8FT4DA5FL1024F
UFA36@G16@OFTFT35@A36@G17@ZFPFPFP4FT1714FZ219D
 
[Main routine]
         T  100 K [define load address]
         G      K [set up relative addressing via @ (theta)]
  [Constants]
     [0] P   10 F [maximum index = 20, edit ad lib.]
     [1] P      D [constant 1]
  [Teleprinter characters]
     [2] #      F [set figures mode]
     [3] C      F [colon (in figures mode)]
     [4] X      F [slash (in figures mode)]
     [5] !      F [space]
     [6] @      F [carriage return]
     [7] &      F [line feed]
     [8] K 4096 F [null]
  [Variables]
     [9] P      F [index]
    [10] P      F [a (where term = a/b)]
    [11] P      F [b]
  [Enter with acc = 0]
    [12] O    2 @ [set teleprinter to figures]
         A    1 @ [acc := 1]
         U    9 @ [index := 1]
         U   10 @ [a := 1]
         T   11 @ [b := 1 (and clear acc)]
         E   34 @ [jump to print first term]
  [Loop back here if not yet printed enough terms]
    [18] A      @ [restore index after test]
         A    1 @ [add 1]
         T    9 @ [update index]
  [Calculate next term. New b := a + b - 2(a mod b).
   Code below calculates c := (a mod b) - b, then new b := a - b - 2*c]
         A   10 @ [acc := a]
    [22] S   11 @ [subtract b]
         E   22 @ [if acc >= 0, subtract again]
         T      F [result c < 0, store in 0F]
         A   10 @ [acc := a]
         S   11 @ [subtract b]
         S      F [subtract c]
         S      F [subtract c]
         T      F [new b = a - b - 2*c; store in 0F]
         A   11 @ [acc := old b]
         T   10 @ [copy to a]
         A      F [acc := new b]
         T   11 @ [copy to b]
  [Print index and a/b. Assume acc = 0 here.]
    [34] A    5 @ [space to replace leading 0's]
         T    1 F [pass to print subroutine]
         A    9 @ [acc := index]
         T      F [pass to print subroutine]
    [38] A   38 @ [for return from subroutine]
         G   56 F [call subroutine, clears acc]
         O    3 @ [print colon]
         O    5 @ [print space]
         A    8 @ [null to replace leading 0's]
         T    1 F [pass to print subroutine]
         A10@ TF A46@ G56F O4@ [print a followed by slash]
         A11@ TF A51@ G56F O6@ O7@ [print b followed by CR LF]
  [Test whether enough terms have been printed]
         A    9 @ [acc := index]
         S      @ [subtract maximum index]
         G   18 @ [loop back if acc < 0]
  [Exit]
         O    8 @ [print null to flush teleprinter buffer]
         Z      F [stop]
         E   12 Z [relative address of entry point]
         P      F [enter with acc = 0]
[end]
 