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_N = 20, max_ciclos = 1000000
 
Function Factorial(Byval N As Integer) As Double
    Dim As Double d: d = 1
    If N = 0 Then Factorial = 1: Exit Function
    While (N > 1)
        d *= N
        N -= 1
    Wend
    Factorial = d
End Function
 
Function Analytical(N As Integer) As Double
    Dim As Double i, sum = 0
    For i = 1 To N
        sum += Factorial(N) / N^i / Factorial(N-i)
    Next i
    Return sum
End Function
 
Function Average(N As Integer, ciclos As Double) As Double
    Dim As Integer i, x, bits, sum = 0
    For i = 0 To ciclos - 1
        x = 1 : bits = 0
        While (bits And x) = 0
            sum += 1
            bits Or= x
            x = 1 Shl (Rnd * (N - 1))
        Wend
    Next i
    Return sum / ciclos
End Function
 
Randomize Timer
Print " N    promedio   analitico    (error)"
Print "--- ---------- ----------- ----------"
For N As Integer = 1 To max_N
    Dim As Double avg = Average(N, max_ciclos)
    Dim As Double ana = Analytical(N)
    Dim As Double diff = abs(avg-ana) / ana * 100
    Print Using " ## #####.###0  #####.###0    ###.#0%"; N; avg; ana; diff
Next N
Sleep
 