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


BEGIN
 
   REAL PROCEDURE FACTORIAL(N); INTEGER N;
   BEGIN
      REAL RESULT;
      INTEGER I;
      RESULT := 1.0;
      FOR I := 2 STEP 1 UNTIL N DO
         RESULT := RESULT * I;
      FACTORIAL := RESULT;
   END FACTORIAL;
 
   REAL PROCEDURE ANALYTICAL (N); INTEGER N;
   BEGIN
      REAL SUM, RN;
      INTEGER I;
      RN := N;
      FOR I := 1 STEP 1 UNTIL N DO
      BEGIN
         SUM := SUM + FACTORIAL(N) / FACTORIAL(N - I) / RN ** I;
      END;
      ANALYTICAL := SUM;
   END ANALYTICAL;
 
   REAL PROCEDURE EXPERIMENTAL(N); INTEGER N;
   BEGIN
      INTEGER NUM;
      INTEGER COUNT;
      INTEGER RUN;
      FOR RUN := 1 STEP 1 UNTIL TESTS DO
      BEGIN
         BOOLEAN ARRAY BITS(1:N);
         INTEGER I;
         FOR I := 1 STEP 1 UNTIL N DO
         BEGIN
            NUM := RANDINT(1,N,SEED);
            IF BITS(NUM) THEN GOTO L;
            BITS(NUM) := TRUE;
            COUNT := COUNT + 1;
         END FOR I;
      L:
      END FOR RUN;
      EXPERIMENTAL := COUNT / TESTS;
   END EXPERIMENTAL;
 
   INTEGER SEED, TESTS;
   SEED := ININT;
   TESTS := 1000000;
   BEGIN
      REAL A, E, ERR;
      INTEGER I;
      OUTTEXT(" N  AVG    CALC   %DIFF"); OUTIMAGE;
      FOR I := 1 STEP 1 UNTIL 20 DO
      BEGIN
         A := ANALYTICAL(I);
         E := EXPERIMENTAL(I);
         ERR := (ABS(E-A)/A)*100.0;
         OUTINT(I, 2);
         OUTFIX(E, 4, 7);
         OUTFIX(A, 4, 10);
         OUTFIX(ERR, 4, 10);
         OUTIMAGE;
      END FOR I;
   END;
END