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


extern crate rand;
 
use rand::{ThreadRng, thread_rng};
use rand::distributions::{IndependentSample, Range};
use std::collections::HashSet;
use std::env;
use std::process;
 
fn help() {
    println!("usage: average_loop_length <max_N> <trials>");
}
 
fn main() {
    let args: Vec<String> = env::args().collect();
    let mut max_n: u32 = 20;
    let mut trials: u32 = 1000;
 
    match args.len() {
        1 => {}
        3 => {
            max_n = args[1].parse::<u32>().unwrap();
            trials = args[2].parse::<u32>().unwrap();
        }
        _ => {
            help();
            process::exit(0);
        }
    }
 
    let mut rng = thread_rng();
 
    println!(" N    average    analytical    (error)");
    println!("===  =========  ============  =========");
    for n in 1..(max_n + 1) {
        let the_analytical = analytical(n);
        let the_empirical = empirical(n, trials, &mut rng);
        println!(" {:>2}     {:3.4}        {:3.4}  ( {:>+1.2}%)",
                 n,
                 the_empirical,
                 the_analytical,
                 100f64 * (the_empirical / the_analytical - 1f64));
    }
}
 
fn factorial(n: u32) -> f64 {
    (1..n + 1).fold(1f64, |p, n| p * n as f64)
}
 
fn analytical(n: u32) -> f64 {
    let sum: f64 = (1..(n + 1))
                       .map(|i| factorial(n) / (n as f64).powi(i as i32) / factorial(n - i))
                       .fold(0f64, |a, v| a + v);
    sum
}
 
fn empirical(n: u32, trials: u32, rng: &mut ThreadRng) -> f64 {
    let sum: f64 = (0..trials)
                       .map(|_t| {
                           let mut item = 1u32;
                           let mut seen = HashSet::new();
                           let range = Range::new(1u32, n + 1);
 
                           for step in 0..n {
                               if seen.contains(&item) {
                                   return step as f64;
                               }
                               seen.insert(item);
                               item = range.ind_sample(rng);
                           }
                           n as f64
                       })
                       .fold(0f64, |a, v| a + v);
    sum / trials as f64
}
 
 
 