Let f be a uniformly-randomly chosen mapping from the numbers 1..N to the numbers 1..N (note: not necessarily a permutation of 1..N; the mapping could produce a number in more than one way or not at all). At some point, the sequence 1, f(1), f(f(1))... will contain a repetition, a number that occurring for the second time in the sequence.



Task

Write a program or a script that estimates, for each N, the average length until the first such repetition.

Also calculate this expected length using an analytical formula, and optionally compare the simulated result with the theoretical one.



This problem comes from the end of Donald Knuth's Christmas tree lecture 2011.

Example of expected output:

 N    average    analytical    (error)
===  =========  ============  =========
  1     1.0000        1.0000  (  0.00%)
  2     1.4992        1.5000  (  0.05%)
  3     1.8784        1.8889  (  0.56%)
  4     2.2316        2.2188  (  0.58%)
  5     2.4982        2.5104  (  0.49%)
  6     2.7897        2.7747  (  0.54%)
  7     3.0153        3.0181  (  0.09%)
  8     3.2429        3.2450  (  0.07%)
  9     3.4536        3.4583  (  0.14%)
 10     3.6649        3.6602  (  0.13%)
 11     3.8091        3.8524  (  1.12%)
 12     3.9986        4.0361  (  0.93%)
 13     4.2074        4.2123  (  0.12%)
 14     4.3711        4.3820  (  0.25%)
 15     4.5275        4.5458  (  0.40%)
 16     4.6755        4.7043  (  0.61%)
 17     4.8877        4.8579  (  0.61%)
 18     4.9951        5.0071  (  0.24%)
 19     5.1312        5.1522  (  0.41%)
 20     5.2699        5.2936  (  0.45%)


/*REXX program computes the average loop length mapping a random field 1···N ───► 1···N */
parse arg runs tests seed .                      /*obtain optional arguments from the CL*/
if  runs =='' |  runs ==","  then runs =      40 /*Not specified?  Then use the default.*/
if tests =='' | tests ==","  then tests= 1000000 /* "      "         "   "   "     "    */
if datatype(seed, 'W')  then call random ,, seed /*Is integer?   For RAND repeatability.*/
!.=0;          !.0=1                             /*used for  factorial (!)  memoization.*/
numeric digits 100000                            /*be able to calculate 25k! if need be.*/
numeric digits max(9, length( !(runs) )   )      /*set the NUMERIC DIGITS for  !(runs). */
say right(     runs, 24)      'runs'             /*display number of runs   we're using.*/
say right(    tests, 24)      'tests'            /*   "       "    " tests    "     "   */
say right( digits(), 24)      'digits'           /*   "       "    " digits   "     "   */
say
say "        N    average     exact     % error "     /* ◄─── title, header ►────────┐  */
hdr="       ═══  ═════════  ═════════  ═════════";       pad=left('',3)  /* ◄────────┘  */
say hdr
         do #=1  for runs;   av=fmtD( exact(#) ) /*use four digits past decimal point.  */
                             xa=fmtD( exper(#) ) /* "    "    "      "     "      "     */
         say right(#,9)  pad xa pad av pad fmtD( abs(xa-av) * 100 / av)   /*show values.*/
         end   /*#*/
say hdr                                          /*display the final header (some bars).*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
!:     procedure expose !.;  parse arg z;                      if !.z\==0  then return !.z
       !=1;       do j=2  for z -1;  !=!*j;  !.j=!;  end; /*compute factorial*/   return !
/*──────────────────────────────────────────────────────────────────────────────────────*/
exact: parse arg x;  s=0;     do j=1  for x;  s=s + !(x) / !(x-j) / x**j;  end;   return s
/*──────────────────────────────────────────────────────────────────────────────────────*/
exper: parse arg n;  k=0;     do tests;   $.=0                      /*do it TESTS times.*/
                                 do n;    r=random(1, n);      if $.r  then leave
                                 $.r=1;   k=k + 1                   /*bump the counter. */
                                 end   /*n*/
                              end      /*tests*/
       return k/tests
/*──────────────────────────────────────────────────────────────────────────────────────*/
fmtD:  parse arg y,d;     d=word(d 4, 1);    y=format(y, , d);     parse var  y  w  '.'  f
       if f=0  then return  w || left('', d +1);                                  return y