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


open System
 
let gamma z = 
    let lanczosCoefficients = [76.18009172947146;-86.50532032941677;24.01409824083091;-1.231739572450155;0.1208650973866179e-2;-0.5395239384953e-5]
    let rec sumCoefficients acc i coefficients =
        match coefficients with
        | []   -> acc
        | h::t -> sumCoefficients (acc + (h/i)) (i+1.0) t
    let gamma = 5.0
    let x = z - 1.0
    Math.Pow(x + gamma + 0.5, x + 0.5) * Math.Exp( -(x + gamma + 0.5) ) * Math.Sqrt( 2.0 * Math.PI ) * sumCoefficients 1.000000000190015 (x + 1.0) lanczosCoefficients
 
let factorial n = gamma ((float n) + 1.)
 
let expected n =
    seq {for i in 1 .. n do yield (factorial n) / System.Math.Pow((float n), (float i)) / (factorial (n - i)) }
    |> Seq.sum
 
let r = System.Random()
 
let trial n =
    let count = ref 0
    let x = ref 1
    let bits = ref 0
    while (!bits &&& !x) = 0 do
        count := !count + 1
        bits := !bits ||| !x
        x := 1 <<< r.Next(n)
    !count
 
 
let tested n times = (float (Seq.sum (seq { for i in 1 .. times do yield (trial n) }))) / (float times)
 
let results = seq {
    for n in 1 .. 20 do
        let avg = tested n 1000000
        let theory = expected n
        yield n, avg, theory
    }
 
 
[<EntryPoint>]
let main argv = 
    printfn " N     average   analytical   (error)"
    printfn "------------------------------------"
    results
    |> Seq.iter (fun (n, avg, theory) ->
        printfn "%2i    %2.6f    %2.6f    %+2.3f%%" n avg theory ((avg / theory - 1.) * 100.))
    0
 