Average loop length: Difference between revisions

Content added Content deleted

Inline

Revision as of 19:43, 3 January 2013

Let f be a randomly chosen mapping from the numbers 1..N to the numbers 1..N. 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.

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

Perl 6

Runs on Rakudo Warszawa (2012.12).

<lang perl6>constant MAX_N = 20; constant TRIALS = 100;

for 1 .. MAX_N -> $N {

   my @lengths = loop-length(random-mapping($N)) xx TRIALS;
   my $average = ([+] @lengths) / @lengths;
   my $analytical = analytical($N);
   my $percent_error = abs($analytical - $average) / $analytical * 100;

   FIRST {
       .say for
       " N    average    analytical    (error)",
       "===  =========  ============  =========";
   }

   printf "%3d  %9.4f  %12.4f    (%4.2f%%)\n",
       $N, $average, $analytical, $percent_error ;

}

sub random-mapping($size) {

   return hash .list Z=> .roll($size) given ^$size;

}

sub loop-length(%mapping) {

   my ($steps, %seen);
   loop (
       my $current = 0;
       !%seen{$current}++;
       $current = %mapping{$current}
   ) { $steps++ }
   return $steps;

}

sub analytical($N) {

   sub postfix:<!>($n) { [*] 1..$n }
   return [+] (1..$N).map: -> $k { $N! * $k * $k / $N ** ($k + 1) / ($N - $k)! };

}</lang>

Example:

$ perl6 loop-lengths
 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%)

@@ Line 37: / Line 37: @@
 Runs on Rakudo Warszawa (2012.12).
-<lang perl6>my $MAX_N  = 20;
+<lang perl6>constant MAX_N  = 20;
-my $TRIALS = 10_000;
+constant TRIALS = 100;
-for 1 .. $MAX_N -> $N {
+for 1 .. MAX_N -> $N {
-    my @lengths = loop-length(random-mapping($N)) xx $TRIALS;
+    my @lengths = loop-length(random-mapping($N)) xx TRIALS;
     my $average = ([+] @lengths) / @lengths;
     my $analytical = analytical($N);
     my $percent_error = abs($analytical - $average) / $analytical * 100;
+    FIRST {
-    FIRST say " N    average    analytical    (error)";
+        .say for
-    FIRST say "===  =========  ============  =========";
+        " N    average    analytical    (error)",
+        "===  =========  ============  =========";
-    my @columns =
+    }
-        $N.fmt('%3d'),
-        $average.fmt('%9.4f'),
-        $analytical.fmt('%12.4f'),
+    printf "%3d  %9.4f  %12.4f    (%4.2f%%)\n",
+        $N, $average, $analytical, $percent_error ;
-        $percent_error.fmt('  (%.2f%%)'),
-    ;
-    say join '  ', @columns;
 }
 sub random-mapping($size) {
+    return hash .list Z=> .roll($size) given ^$size;
-    # Not just a permutation, but a mapping from all 1..N to any 1..N
-    # In other words, .roll, not .pick
-    return hash map { $_ => (1..$size).roll }, 1..$size;
 }
 sub loop-length(%mapping) {
-    my %seen;
+    my ($steps, %seen);
-    my $current = 1;
+    loop (
-    my $steps = 0;
+        my $current = 0;
-    while !%seen{$current}++ {
+        !%seen{$current}++;
-        $current = %mapping{$current};
+        $current = %mapping{$current}
-        $steps++;
+    ) { $steps++ }
-    }
     return $steps;
 }