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%)


 
Const MAX = 20
Const ITER = 100000
 
Function expected(n) 
  Dim sum
  ni=n
  For i = 1 To n
    sum = sum + fact(n) / ni / fact(n-i)
    ni=ni*n
  Next
  expected = sum
End Function
 
Function test(n )
  Dim coun,x,bits 
  For i = 1 To ITER
    x = 1
    bits = 0
    Do While Not bits And x
      count = count + 1
      bits = bits Or x
      x =shift(Int(n * Rnd()))
    Loop
  Next
  test = count / ITER
End Function
 
'VBScript formats numbers but does'nt align them!
function rf(v,n,s) rf=right(string(n,s)& v,n):end function  
 
'some precalculations to speed things up...    
dim fact(20),shift(20)
fact(0)=1:shift(0)=1
for i=1 to 20
  fact(i)=i*fact(i-1)
  shift(i)=2*shift(i-1)
next
 
Dim n 
Wscript.echo "For " & ITER &" iterations" 
Wscript.Echo " n     avg.     exp.   (error%)"
Wscript.Echo "==   ======   ======  =========="
For n = 1 To MAX
  av = test(n)
  ex = expected(n)
  Wscript.Echo rf(n,2," ")& "  "& rf(formatnumber(av, 4),7," ") & "   "& _
  rf(formatnumber(ex,4),6," ")& "  ("& rf(Formatpercent(1 - av / ex,4),8," ") & ")"
Next
 