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


/*REXX pgm finds the Nth value of the  Calkin─Wilf  sequence (which will be a fraction),*/
/*────────────────────── or finds which sequence number contains a specified fraction). */
numeric digits 2000                              /*be able to handle ginormic integers. */
parse arg LO HI te .                             /*obtain optional arguments from the CL*/
if LO=='' | LO==","   then LO=  1                /*Not specified?  Then use the default.*/
if HI=='' | HI==","   then HI= 20                /* "      "         "   "   "     "    */
if te=='' | te==","   then te= '/'               /* "      "         "   "   "     "    */
if datatype(LO, 'W')  then call CW_terms         /*Is LO numeric?  Then show some terms.*/
if pos('/', te)>0     then call CW_frac  te      /*Does TE have a / ?   Then find term #*/
exit 0
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg ?;  do jc=length(?)-3  to 1  by -3; ?=insert(',', ?, jc); end;  return ?
th:     parse arg th; return word('th st nd rd', 1+(th//10) *(th//100%10\==1) *(th//10<4))
/*──────────────────────────────────────────────────────────────────────────────────────*/
CW_frac:   procedure; parse arg p '/' q .;       say
           if q==''  then do;  p= 83116;         q= 51639;  end
           n= rle2dec( frac2cf(p q) );                    @CWS= 'the Calkin─Wilf sequence'
           say 'for '  p"/"q',  the element number for'   @CWS    "is: "    commas(n)th(n)
           if length(n)<10  then return
           say;  say 'The above number has '     commas(length(n))      " decimal digits."
           return
/*──────────────────────────────────────────────────────────────────────────────────────*/
CW_term:   procedure;  parse arg z;                 dd= 1;               nn= 0
                                       do z
                                       parse value  dd  dd*(2*(nn%dd)+1)-nn   with  nn  dd
                                       end   /*z*/
           return nn'/'dd
/*──────────────────────────────────────────────────────────────────────────────────────*/
CW_terms:  $=;        if LO\==0  then  do j=LO  to HI;   $= $  CW_term(j)','
                                       end   /*j*/
           if $==''  then return
           say 'Calkin─Wilf sequence terms for '  commas(LO)  " ──► "  commas(HI)  ' are:'
           say strip( strip($), 'T', ",")
           return
/*──────────────────────────────────────────────────────────────────────────────────────*/
frac2cf:   procedure;  parse arg p q;  if q==''  then return p;          cf= p % q;   m= q
           p= p - cf*q;                n= p;        if p==0  then return cf
                         do k=1  until n==0;        @.k= m % n
                         m= m  -  @.k * n;    parse value  n m   with   m n   /*swap N M*/
                         end   /*k*/
                                              /*for inverse Calkin─Wilf, K must be even.*/
           if k//2  then do;  @.k= @.k - 1;   k= k + 1;    @.k= 1;   end
                         do k=1  for k;       cf= cf @.k;            end  /*k*/
           return cf
/*──────────────────────────────────────────────────────────────────────────────────────*/
rle2dec:   procedure;  parse arg f1 rle;                       obin= copies(1, f1)
                               do until rle=='';               parse var rle f0 f1 rle
                               obin= copies(1, f1)copies(0, f0)obin
                               end   /*until*/
           return x2d( b2x(obin) )            /*RLE2DEC: Run Length Encoding ──► decimal*/