CloudFlare suffered a massive security issue affecting all of its customers, including Rosetta Code. All passwords not changed since February 19th 2017 have been expired, and session cookie longevity will be reduced until late March.--Michael Mol (talk) 05:15, 25 February 2017 (UTC)

Average loop length

From Rosetta Code
Task
Average loop length
You are encouraged to solve this task according to the task description, using any language you may know.

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


Ada[edit]

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Numerics.Discrete_Random;
procedure Avglen is
package IIO is new Ada.Text_IO.Integer_IO (Positive); use IIO;
package LFIO is new Ada.Text_IO.Float_IO (Long_Float); use LFIO;
subtype FactN is Natural range 0..20;
TESTS : constant Natural := 1_000_000;
 
function Factorial (N : FactN) return Long_Float is
Result : Long_Float := 1.0;
begin
for I in 2..N loop Result := Result * Long_Float(I); end loop;
return Result;
end Factorial;
 
function Analytical (N : FactN) return Long_Float is
Sum : Long_Float := 0.0;
begin
for I in 1..N loop
Sum := Sum + Factorial(N) / Factorial(N - I) / Long_Float(N)**I;
end loop;
return Sum;
end Analytical;
 
function Experimental (N : FactN) return Long_Float is
subtype RandInt is Natural range 1..N;
package Random is new Ada.Numerics.Discrete_Random(RandInt);
seed : Random.Generator;
Num : RandInt;
count : Natural := 0;
bits : array(RandInt'Range) of Boolean;
begin
Random.Reset(seed);
for run in 1..TESTS loop
bits := (others => false);
for I in RandInt'Range loop
Num := Random.Random(seed); exit when bits(Num);
bits(Num) := True; count := count + 1;
end loop;
end loop;
return Long_Float(count)/Long_Float(TESTS);
end Experimental;
 
A, E, err : Long_Float;
begin
Put_Line(" N avg calc  %diff");
for I in 1..20 loop
A := Analytical(I); E := Experimental(I); err := abs(E-A)/A*100.0;
Put(I, Width=>2); Put(E ,Aft=>4, exp=>0); Put(A, Aft=>4, exp=>0);
Put(err, Fore=>3, Aft=>3, exp=>0); New_line;
end loop;
end Avglen;
Output:
 N  avg    calc   %diff
 1 1.0000 1.0000  0.000
 2 1.5000 1.5000  0.003
 3 1.8886 1.8889  0.015
 4 2.2180 2.2188  0.033
 5 2.5104 2.5104  0.000
 6 2.7745 2.7747  0.006
 7 3.0191 3.0181  0.033
 8 3.2433 3.2450  0.052
 9 3.4583 3.4583  0.001
10 3.6597 3.6602  0.015
11 3.8524 3.8524  0.001
12 4.0352 4.0361  0.022
13 4.2147 4.2123  0.055
14 4.3853 4.3820  0.075
15 4.5453 4.5458  0.011
16 4.7055 4.7043  0.027
17 4.8592 4.8579  0.028
18 5.0062 5.0071  0.017
19 5.1535 5.1522  0.025
20 5.2955 5.2936  0.035

BBC BASIC[edit]

      @% = &2040A
MAX_N = 20
TIMES = 1000000
 
FOR n = 1 TO MAX_N
avg = FNtest(n, TIMES)
theory = FNanalytical(n)
diff = (avg / theory - 1) * 100
PRINT STR$(n), avg, theory, diff "%"
NEXT
END
 
DEF FNanalytical(n)
LOCAL i, s
FOR i = 1 TO n
s += FNfactorial(n) / n^i / FNfactorial(n-i)
NEXT
= s
 
DEF FNtest(n, times)
LOCAL i, b, c, x
FOR i = 1 TO times
x = 1 : b = 0
WHILE (b AND x) = 0
c += 1
b OR= x
x = 1 << (RND(n) - 1)
ENDWHILE
NEXT
= c / times
 
DEF FNfactorial(n)
IF n=1 OR n=0 THEN =1 ELSE = n * FNfactorial(n-1)
Output:
1             1.0000    1.0000    0.0000%
2             1.4995    1.5000   -0.0366%
3             1.8879    1.8889   -0.0509%
4             2.2193    2.2188    0.0240%
5             2.5105    2.5104    0.0057%
6             2.7755    2.7747    0.0293%
7             3.0199    3.0181    0.0573%
8             3.2396    3.2450   -0.1664%
9             3.4562    3.4583   -0.0609%
10            3.6578    3.6602   -0.0659%
11            3.8523    3.8524   -0.0025%
12            4.0336    4.0361   -0.0602%
13            4.2139    4.2123    0.0366%
14            4.3816    4.3820   -0.0105%
15            4.5432    4.5458   -0.0570%
16            4.7108    4.7043    0.1386%
17            4.8578    4.8579   -0.0018%
18            5.0063    5.0071   -0.0144%
19            5.1564    5.1522    0.0814%
20            5.2945    5.2936    0.0166%

C[edit]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
 
#define MAX_N 20
#define TIMES 1000000
 
double factorial(int n) {
double f = 1;
int i;
for (i = 1; i <= n; i++) f *= i;
return f;
}
 
double expected(int n) {
double sum = 0;
int i;
for (i = 1; i <= n; i++)
sum += factorial(n) / pow(n, i) / factorial(n - i);
return sum;
}
 
int randint(int n) {
int r, rmax = RAND_MAX / n * n;
while ((r = rand()) >= rmax);
return r / (RAND_MAX / n);
}
 
int test(int n, int times) {
int i, count = 0;
for (i = 0; i < times; i++) {
int x = 1, bits = 0;
while (!(bits & x)) {
count++;
bits |= x;
x = 1 << randint(n);
}
}
return count;
}
 
int main(void) {
srand(time(0));
puts(" n\tavg\texp.\tdiff\n-------------------------------");
 
int n;
for (n = 1; n <= MAX_N; n++) {
int cnt = test(n, TIMES);
double avg = (double)cnt / TIMES;
double theory = expected(n);
double diff = (avg / theory - 1) * 100;
printf("%2d %8.4f %8.4f %6.3f%%\n", n, avg, theory, diff);
}
return 0;
}
Output:
 n      avg     exp.    diff
-------------------------------
 1   1.0000   1.0000  0.000%
 2   1.4998   1.5000 -0.015%
 3   1.8879   1.8889 -0.051%
 4   2.2181   2.2188 -0.029%
 5   2.5107   2.5104  0.012%
 6   2.7741   2.7747 -0.021%
 7   3.0168   3.0181 -0.044%
 8   3.2455   3.2450  0.014%
 9   3.4591   3.4583  0.023%
10   3.6596   3.6602 -0.017%
11   3.8519   3.8524 -0.013%
12   4.0384   4.0361  0.059%
13   4.2106   4.2123 -0.042%
14   4.3840   4.3820  0.044%
15   4.5449   4.5458 -0.020%
16   4.7058   4.7043  0.033%
17   4.8549   4.8579 -0.060%
18   5.0084   5.0071  0.026%
19   5.1479   5.1522 -0.084%
20   5.2957   5.2936  0.040%

Clojure[edit]

Translation of: Python
(ns cyclelengths
(:gen-class))
 
(defn factorial [n]
" n! "
(apply *' (range 1 (inc n)))) ; Use *' (vs. *) to allow arbitrary length arithmetic
 
(defn pow [n i]
" n^i"
(apply *' (repeat i n)))
 
(defn analytical [n]
" Analytical Computation "
(->>(range 1 (inc n))
(map #(/ (factorial n) (pow n %) (factorial (- n %)))) ;calc n %))
(reduce + 0)))
 
;; Number of random times to test each n
(def TIMES 1000000)
 
(defn single-test-cycle-length [n]
" Single random test of cycle length "
(loop [count 0
bits 0
x 1]
(if (zero? (bit-and x bits))
(recur (inc count) (bit-or bits x) (bit-shift-left 1 (rand-int n)))
count)))
 
(defn avg-cycle-length [n times]
" Average results of single tests of cycle lengths "
(/
(reduce +
(for [i (range times)]
(single-test-cycle-length n)))
times))
 
;; Show Results
(println "\tAvg\t\tExp\t\tDiff")
(doseq [q (range 1 21)
:let [anal (double (analytical q))
avg (double (avg-cycle-length q TIMES))
diff (Math/abs (* 100 (- 1 (/ avg anal))))]]
(println (format "%3d\t%.4f\t%.4f\t%.2f%%" q avg anal diff)))
 
Output:
	Avg	Exp	Diff
  1	1.0000	1.0000	0.00%
  2	1.4995	1.5000	0.03%
  3	1.8899	1.8889	0.05%
  4	2.2178	2.2188	0.04%
  5	2.5118	2.5104	0.06%
  6	2.7773	2.7747	0.09%
  7	3.0177	3.0181	0.02%
  8	3.2448	3.2450	0.01%
  9	3.4587	3.4583	0.01%
 10	3.6594	3.6602	0.02%
 11	3.8553	3.8524	0.08%
 12	4.0335	4.0361	0.06%
 13	4.2113	4.2123	0.03%
 14	4.3823	4.3820	0.01%
 15	4.5491	4.5458	0.07%
 16	4.7035	4.7043	0.02%
 17	4.8580	4.8579	0.00%
 18	5.0050	5.0071	0.04%
 19	5.1543	5.1522	0.04%
 20	5.2956	5.2936	0.04%

D[edit]

Translation of: Perl 6
import std.stdio, std.random, std.math, std.algorithm, std.range, std.format;
 
real analytical(in int n) pure nothrow @safe [email protected]*/ {
enum aux = (int k) => reduce!q{a * b}(1.0L, iota(n - k + 1, n + 1));
return iota(1, n + 1)
.map!(k => (aux(k) * k ^^ 2) / (real(n) ^^ (k + 1)))
.sum;
}
 
size_t loopLength(size_t maxN)(in int size, ref Xorshift rng) {
__gshared static bool[maxN + 1] seen;
seen[0 .. size + 1] = false;
int current = 1;
size_t steps = 0;
while (!seen[current]) {
seen[current] = true;
current = uniform(1, size + 1, rng);
steps++;
}
return steps;
}
 
void main() {
enum maxN = 40;
enum nTrials = 300_000;
auto rng = Xorshift(unpredictableSeed);
writeln(" n average analytical (error)");
writeln("=== ========= ============ ==========");
 
foreach (immutable n; 1 .. maxN + 1) {
long total = 0;
foreach (immutable _; 0 .. nTrials)
total += loopLength!maxN(n, rng);
immutable average = total / real(nTrials);
immutable an = n.analytical;
immutable percentError = abs(an - average) / an * 100;
immutable errorS = format("%2.4f", percentError);
writefln("%3d  %9.5f  %12.5f (%7s%%)",
n, average, an, errorS);
}
}
Output:
 n    average    analytical     (error)
===  =========  ============  ==========
  1    1.00000       1.00000  ( 0.0000%)
  2    1.50017       1.50000  ( 0.0111%)
  3    1.88932       1.88889  ( 0.0226%)
  4    2.21795       2.21875  ( 0.0362%)
  5    2.51159       2.51040  ( 0.0474%)
  6    2.77373       2.77469  ( 0.0345%)
  7    3.01894       3.01814  ( 0.0264%)
  8    3.24734       3.24502  ( 0.0716%)
  9    3.45876       3.45832  ( 0.0127%)
 10    3.66595       3.66022  ( 0.1567%)
 11    3.85000       3.85237  ( 0.0616%)
 12    4.03532       4.03607  ( 0.0187%)
 13    4.20879       4.21235  ( 0.0843%)
 14    4.37664       4.38203  ( 0.1230%)
 15    4.54986       4.54581  ( 0.0892%)
 16    4.70431       4.70426  ( 0.0010%)
 17    4.85640       4.85787  ( 0.0302%)
 18    5.01359       5.00706  ( 0.1303%)
 19    5.15487       5.15220  ( 0.0519%)
 20    5.29486       5.29358  ( 0.0241%)
 21    5.43276       5.43150  ( 0.0231%)
 22    5.56570       5.56620  ( 0.0088%)
 23    5.70611       5.69788  ( 0.1443%)
 24    5.82618       5.82675  ( 0.0098%)
 25    5.94846       5.95298  ( 0.0759%)
 26    6.07440       6.07672  ( 0.0381%)
 27    6.20717       6.19811  ( 0.1461%)
 28    6.31546       6.31729  ( 0.0290%)
 29    6.44201       6.43437  ( 0.1187%)
 30    6.54592       6.54946  ( 0.0540%)
 31    6.65818       6.66265  ( 0.0671%)
 32    6.77215       6.77405  ( 0.0279%)
 33    6.88381       6.88372  ( 0.0013%)
 34    6.99790       6.99175  ( 0.0880%)
 35    7.10990       7.09820  ( 0.1648%)
 36    7.20391       7.20316  ( 0.0104%)
 37    7.30085       7.30667  ( 0.0796%)
 38    7.40366       7.40880  ( 0.0693%)
 39    7.51864       7.50959  ( 0.1204%)
 40    7.60255       7.60911  ( 0.0863%)

EchoLisp[edit]

 
(lib 'math) ;; Σ aka (sigma f(n) nfrom nto)
 
(define (f-count N (times 100000))
(define count 0)
(for ((i times))
 
;; new random f mapping from 0..N-1 to 0..N-1
;; (f n) is NOT (random N)
;; because each call (f n) must return the same value
 
(define f (build-vector N (lambda(i) (random N))))
 
(define hits (make-vector N))
(define n 0)
(while (zero? [hits n])
(++ count)
(vector+= hits n 1)
(set! n [f n])))
(// count times))
 
(define (f-anal N)
(Σ (lambda(i) (// (! N) (! (- N i)) (^ N i))) 1 N))
 
(decimals 5)
(define (f-print (maxN 21))
(for ((N (in-range 1 maxN)))
(define fc (f-count N))
(define fa (f-anal N))
(printf "%3d %10d %10d  %10.2d %%" N fc fa (// (abs (- fa fc)) fc 0.01))))
 
Output:
(f-print)
  1          1          1          0 %
  2    1.49908        1.5       0.06 %
  3    1.89059    1.88889       0.09 %
  4    2.21709    2.21875       0.07 %
  5    2.50629     2.5104       0.16 %
  6    2.77027    2.77469       0.16 %
  7    3.01739    3.01814       0.02 %
  8    3.23934    3.24502       0.18 %
  9    3.45862    3.45832       0.01 %
 10    3.65959    3.66022       0.02 %
 11    3.85897    3.85237       0.17 %
 12    4.04188    4.03607       0.14 %
 13    4.21226    4.21235          0 %
 14    4.38021    4.38203       0.04 %
 15    4.54158    4.54581       0.09 %
 16    4.70633    4.70426       0.04 %
 17    4.86109    4.85787       0.07 %
 18    4.99903    5.00706       0.16 %
 19    5.15873     5.1522       0.13 %
 20    5.30243    5.29358       0.17 %
 


Elixir[edit]

Translation of: Ruby
Works with: Elixir version 1.1+
defmodule RC do
def factorial(0), do: 1
def factorial(n), do: Enum.reduce(1..n, 1, &(&1 * &2))
 
def loop_length(n), do: loop_length(n, MapSet.new)
 
defp loop_length(n, set) do
r = :rand.uniform(n)
if r in set, do: MapSet.size(set), else: loop_length(n, MapSet.put(set, r))
end
 
def task(runs) do
IO.puts " N average analytical (error) "
IO.puts "=== ========= ========== ========="
Enum.each(1..20, fn n ->
avg = Enum.reduce(1..runs, 0, fn _,sum -> sum + loop_length(n) end) / runs
analytical = Enum.reduce(1..n, 0, fn i,sum ->
sum + (factorial(n) / :math.pow(n, i) / factorial(n-i))
end)
 :io.format "~3w ~9.4f ~9.4f (~6.2f%)~n", [n, avg, analytical, abs(avg/analytical - 1)*100]
end)
end
end
 
runs = 1_000_000
RC.task(runs)
Output:
 N    average   analytical   (error)
===  =========  ==========  =========
  1     1.0000      1.0000  (  0.00%)
  2     1.5001      1.5000  (  0.00%)
  3     1.8892      1.8889  (  0.02%)
  4     2.2189      2.2188  (  0.01%)
  5     2.5113      2.5104  (  0.04%)
  6     2.7749      2.7747  (  0.01%)
  7     3.0185      3.0181  (  0.01%)
  8     3.2456      3.2450  (  0.02%)
  9     3.4612      3.4583  (  0.08%)
 10     3.6573      3.6602  (  0.08%)
 11     3.8524      3.8524  (  0.00%)
 12     4.0357      4.0361  (  0.01%)
 13     4.2102      4.2123  (  0.05%)
 14     4.3813      4.3820  (  0.02%)
 15     4.5422      4.5458  (  0.08%)
 16     4.7057      4.7043  (  0.03%)
 17     4.8581      4.8579  (  0.01%)
 18     5.0045      5.0071  (  0.05%)
 19     5.1533      5.1522  (  0.02%)
 20     5.2951      5.2936  (  0.03%)

F#[edit]

Translation of: Scala

But uses the Gamma function instead of factorials.

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
 
Output:
 N     average   analytical   (error)
------------------------------------
 1    1.000000    1.000000    +0.000%
 2    1.498934    1.500000    -0.071%
 3    1.889318    1.888889    +0.023%
 4    2.219397    2.218750    +0.029%
 5    2.510618    2.510400    +0.009%
 6    2.771914    2.774691    -0.100%
 7    3.014726    3.018139    -0.113%
 8    3.245022    3.245018    +0.000%
 9    3.457096    3.458316    -0.035%
10    3.660337    3.660216    +0.003%
11    3.849770    3.852372    -0.068%
12    4.038977    4.036074    +0.072%
13    4.213248    4.212348    +0.021%
14    4.380451    4.382029    -0.036%
15    4.541868    4.545807    -0.087%
16    4.704117    4.704258    -0.003%
17    4.858934    4.857871    +0.022%
18    5.004236    5.007063    -0.056%
19    5.154166    5.152196    +0.038%
20    5.298119    5.293585    +0.086%

Go[edit]

package main
 
import (
"fmt"
"math"
"math/rand"
)
 
const nmax = 20
 
func main() {
fmt.Println(" N average analytical (error)")
fmt.Println("=== ========= ============ =========")
for n := 1; n <= nmax; n++ {
a := avg(n)
b := ana(n)
fmt.Printf("%3d  %9.4f  %12.4f (%6.2f%%)\n",
n, a, b, math.Abs(a-b)/b*100)
}
}
 
func avg(n int) float64 {
const tests = 1e4
sum := 0
for t := 0; t < tests; t++ {
var v [nmax]bool
for x := 0; !v[x]; x = rand.Intn(n) {
v[x] = true
sum++
}
}
return float64(sum) / tests
}
 
func ana(n int) float64 {
nn := float64(n)
term := 1.
sum := 1.
for i := nn - 1; i >= 1; i-- {
term *= i / nn
sum += term
}
return sum
}
Output:
 N    average    analytical    (error)
===  =========  ============  =========
  1     1.0000        1.0000  (  0.00%)
  2     1.5007        1.5000  (  0.05%)
  3     1.8959        1.8889  (  0.37%)
  4     2.2138        2.2188  (  0.22%)
  5     2.5013        2.5104  (  0.36%)
  6     2.7940        2.7747  (  0.70%)
  7     3.0197        3.0181  (  0.05%)
  8     3.2715        3.2450  (  0.82%)
  9     3.4147        3.4583  (  1.26%)
 10     3.6758        3.6602  (  0.43%)
 11     3.8672        3.8524  (  0.38%)
 12     4.0309        4.0361  (  0.13%)
 13     4.2153        4.2123  (  0.07%)
 14     4.3380        4.3820  (  1.00%)
 15     4.5030        4.5458  (  0.94%)
 16     4.7563        4.7043  (  1.11%)
 17     4.8616        4.8579  (  0.08%)
 18     4.9933        5.0071  (  0.27%)
 19     5.1534        5.1522  (  0.02%)
 20     5.3031        5.2936  (  0.18%)

Haskell[edit]

import System.Random
import qualified Data.Set as S
import Text.Printf
 
findRep :: (Random a, Integral a, RandomGen b) => a -> b -> (a, b)
findRep n gen = findRep' (S.singleton 1) 1 gen
where
findRep'
seen len gen'
| S.member fx seen = (len, gen'
')
| otherwise = findRep'
(S.insert fx seen) (len + 1) gen''
where
(fx, gen'') = randomR (1, n) gen'
 
statistical :: (Integral a, Random b, Integral b, RandomGen c, Fractional d) =>
a -> b -> c -> (d, c)
statistical samples size gen =
let (total, gen'
) = sar samples gen 0
in ((fromIntegral total) / (fromIntegral samples), gen')
where
sar 0 gen'
acc = (acc, gen')
sar samples'
gen' acc =
let (len, gen'
') = findRep size gen'
in sar (samples' - 1) gen'' (acc + len)
 
factorial :: (Integral a) => a -> a
factorial n = foldl (*) 1 [1..n]
 
analytical :: (Integral a, Fractional b) => a -> b
analytical n = sum [fromIntegral num /
fromIntegral (factorial (n - i)) /
fromIntegral (n ^ i) |
i <- [1..n]]
where num = factorial n
 
test :: (Integral a, Random b, Integral b, PrintfArg b, RandomGen c) =>
a -> [b] -> c -> IO c
test _ [] gen = return gen
test samples (x:xs) gen = do
let (st, gen'
) = statistical samples x gen
an = analytical x
err = abs (st - an) / st * 100.0
str = printf "%3d  %9.4f  %12.4f (%6.2f%%)\n"
x (st :: Float) (an :: Float) (err :: Float)
putStr str
test samples xs gen'
 
main :: IO ()
main = do
putStrLn " N average analytical (error)"
putStrLn "=== ========= ============ ========="
let samples = 10000 :: Integer
range = [1..20] :: [Integer]
_ <- test samples range $ mkStdGen 0
return ()
 N    average    analytical    (error)
===  =========  ============  =========
  1     1.0000        1.0000  (  0.00%)
  2     1.4941        1.5000  (  0.39%)
  3     1.8895        1.8889  (  0.03%)
  4     2.2246        2.2188  (  0.26%)
  5     2.5158        2.5104  (  0.21%)
  6     2.7875        2.7747  (  0.46%)
  7     3.0425        3.0181  (  0.80%)
  8     3.2157        3.2450  (  0.91%)
  9     3.4534        3.4583  (  0.14%)
 10     3.6561        3.6602  (  0.11%)
 11     3.8357        3.8524  (  0.43%)
 12     4.0291        4.0361  (  0.17%)
 13     4.1819        4.2123  (  0.73%)
 14     4.3469        4.3820  (  0.81%)
 15     4.4942        4.5458  (  1.15%)
 16     4.7093        4.7043  (  0.11%)
 17     4.8288        4.8579  (  0.60%)
 18     5.0021        5.0071  (  0.10%)
 19     5.1980        5.1522  (  0.88%)
 20     5.2961        5.2936  (  0.05%)

J[edit]

First, let's consider an exact, brute force approach.

Since J array indices start at 0, we'll work with 0..N-1 instead of 1..N, dealing with the difference at the boundaries.

We can implement f as {&LIST where LIST is an arbitrary list of N numbers, each picked independently from the range 0..(N-1). We can incrementally build the described sequence using (, f@{:) - here we extend the sequence by applying f to the last element of the sequence. Since we are only concerned with the sequence up to the point of the first repeat, we can select the unique values giving us (~.@, f@{:). This routine stops changing when we reach the desired length, so we can repeatedly apply it forever. For example:

   (~.@, {&0 0@{:)^:_] 0
0
(~.@, {&0 0@{:)^:_] 1
1 0

Once we have the sequence, we can count how many elements are in it.

   0 0 ([: # (] ~.@, {:@] { [)^:_) 1
2

Meanwhile, we can also generate all possible values of 1..N by counting out N^N values and breaking out the result as a base N list of digits.

   (#.inv i.@^~)2
0 0
0 1
1 0
1 1

All that's left is to count the lengths of all possible sequences for all possible distinct instances of f and average the results:

   (+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)1
1
(+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)2
1.5
(+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)3
1.88889
(+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)4
2.21875
(+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)5
2.5104
(+/ % #)@,@((#.inv i.@^~) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)6
2.77469

Meanwhile the analytic solution (derived by reading the Ada implementation) looks like this:

   ana=: +/@(!@[ % !@- * ^) 1+i.
ana"0]1 2 3 4 5 6
1 1.5 1.88889 2.21875 2.5104 2.77469

To get our simulation, we can take the exact approach and replace the part that generates all possible values for f with a random mechanism. Since the task does not specify how long to run the simulation, and to make this change easy, we'll use N*1e4 tests.

   sim=: (+/ % #)@,@((]?@$~1e4,]) ([: # (] ~.@, {:@] { [)^:_)"1 0/ i.)
sim"0]1 2 3 4 5 6
1 1.5034 1.8825 2.22447 2.51298 2.76898

The simulation approach is noticeably slower than the analytic approach, while being less accurate.

Finally, we can generate our desired results:

   (;:'N average analytic error'),:,.each(;ana"0 ([;];-|@%[) sim"0)1+i.20
+--+-------+--------+-----------+
|N |average|analytic|error |
+--+-------+--------+-----------+
| 1| 1| 1 | 0|
| 2| 1.5|1.49955 | 0.0003|
| 3|1.88889| 1.8928 | 0.00207059|
| 4|2.21875|2.23082 | 0.00544225|
| 5| 2.5104|2.52146 | 0.00440567|
| 6|2.77469|2.78147 | 0.00244182|
| 7|3.01814| 3.0101 | 0.00266346|
| 8|3.24502|3.25931 | 0.00440506|
| 9|3.45832|3.45314 | 0.00149532|
|10|3.66022| 3.6708 | 0.00289172|
|11|3.85237|3.84139 | 0.00285049|
|12|4.03607|4.03252 |0.000881304|
|13|4.21235|4.18358 | 0.00682833|
|14|4.38203|4.38791 | 0.00134132|
|15|4.54581|4.54443 |0.000302246|
|16|4.70426|4.71351 | 0.00196721|
|17|4.85787|4.85838 |0.000104089|
|18|5.00706|5.00889 |0.000365752|
|19| 5.1522|5.14785 |0.000843052|
|20|5.29358|5.28587 | 0.00145829|
+--+-------+--------+-----------+

Here, error is the difference between the two values divided by the analytic value.

Java[edit]

This uses a 0-based index (0, 1, ..., n-1) as opposed to the 1-based index (1, 2, ..., n) specified in the question, because it fits better with the native structure of Java.

import java.util.ArrayList;
 
public class AverageLoopLength {
private static final int N = 100000;
//analytical(n) = sum_(i=1)^n (n!/(n-i)!/n**i)
public static float analytical(int n){
float[] factorial = new float[n+1];
float[] powers = new float[n+1];
factorial[0] = powers[0] = 1;
for(int i=1;i<=n;i++){
factorial[i] = factorial[i-1] * i;
powers[i] = powers[i-1] * n;
}
float sum = 0;
//memoized factorial and powers
for(int i=1;i<=n;i++){
sum += factorial[n]/factorial[n-i]/powers[i];
}
return sum;
}
public static float average(int n){
float sum = 0;
for(int a=0;a<N;a++){
int[] random = new int[n];
for(int i=0;i<n;i++){
random[i] = (int)(Math.random()*n);
}
ArrayList<Integer> seen = new ArrayList<>(n);
int current = 0;
int length = 0;
while(true){
length++;
seen.add(current);
current = random[current];
if(seen.contains(current)){
break;
}
}
sum += length;
}
return sum/N;
}
public static void main(String args[]){
System.out.println(" N average analytical (error)\n=== ========= ============ =========");
for(int i=1;i<=20;i++){
float avg = average(i);
float ana = analytical(i);
System.out.println(String.format("%3d  %9.4f  %12.4f (%6.2f%%)",i,avg,ana,((ana-avg)/ana*100)));;
}
}
}
 

Kotlin[edit]

Translation of: Go
// version 1.0.6
 
const val NMAX = 20
const val TESTS = 1000000
val rand = java.util.Random()
 
fun avg(n: Int): Double {
var sum = 0
for (t in 0 until TESTS) {
val v = BooleanArray(NMAX)
var x = 0
while (!v[x]) {
v[x] = true
sum++
x = rand.nextInt(n)
}
}
return sum.toDouble() / TESTS
}
 
fun ana(n: Int): Double {
var nn = n.toDouble()
var term = 1.0
var sum = 1.0
for (i in n - 1 downTo 1) {
term *= i / nn
sum += term
}
return sum
}
 
fun main(args: Array<String>) {
println(" N average analytical (error)")
println("=== ========= ============ =========")
for (n in 1..NMAX) {
val a = avg(n)
val b = ana(n)
println(String.format("%3d  %6.4f  %10.4f (%4.2f%%)", n, a, b, Math.abs(a - b) / b * 100.0))
}
}

Sample output:

Output:
 N    average    analytical    (error)
===  =========  ============  =========
  1   1.0000       1.0000      (0.00%)
  2   1.5004       1.5000      (0.03%)
  3   1.8890       1.8889      (0.00%)
  4   2.2179       2.2188      (0.04%)
  5   2.5108       2.5104      (0.02%)
  6   2.7738       2.7747      (0.03%)
  7   3.0178       3.0181      (0.01%)
  8   3.2482       3.2450      (0.10%)
  9   3.4572       3.4583      (0.03%)
 10   3.6608       3.6602      (0.02%)
 11   3.8545       3.8524      (0.06%)
 12   4.0378       4.0361      (0.04%)
 13   4.2131       4.2123      (0.02%)
 14   4.3795       4.3820      (0.06%)
 15   4.5481       4.5458      (0.05%)
 16   4.7044       4.7043      (0.00%)
 17   4.8610       4.8579      (0.06%)
 18   5.0027       5.0071      (0.09%)
 19   5.1498       5.1522      (0.05%)
 20   5.2941       5.2936      (0.01%)

Liberty BASIC[edit]

Translation of: BBC BASIC
 
MAXN = 20
TIMES = 10000'00
 
't0=time$("ms")
FOR n = 1 TO MAXN
avg = FNtest(n, TIMES)
theory = FNanalytical(n)
diff = (avg / theory - 1) * 100
PRINT n, avg, theory, using("##.####",diff); "%"
NEXT
't1=time$("ms")
'print t1-t0; " ms"
END
 
function FNanalytical(n)
FOR i = 1 TO n
s = s+ FNfactorial(n) / n^i / FNfactorial(n-i)
NEXT
FNanalytical = s
end function
 
function FNtest(n, times)
FOR i = 1 TO times
x = 1 : b = 0
WHILE (b AND x) = 0
c = c + 1
b = b OR x
x = 2^int(n*RND(1))
WEND
NEXT
FNtest = c / times
end function
 
function FNfactorial(n)
IF n=1 OR n=0 THEN FNfactorial=1 ELSE FNfactorial= n * FNfactorial(n-1)
end function
 
Output:
1             1             1              0.0000%
2             1.4759        1.5           -1.6067%
3             1.8868        1.88888889    -0.1106%
4             2.2139        2.21875       -0.2186%
5             2.4784        2.5104        -1.2747%
6             2.7888        2.77469136     0.5085%
7             2.9846        3.0181387     -1.1112%
8             3.2645        3.24501801     0.6004%
9             3.464         3.45831574     0.1644%
10            3.6602        3.66021568    -0.0004%
11            3.8255        3.85237205    -0.6975%
12            4.019         4.03607368    -0.4230%
13            4.2033        4.21234791    -0.2148%
14            4.3985        4.38202942     0.3759%
15            4.5868        4.54580729     0.9018%
16            4.6705        4.70425825    -0.7176%
17            4.8807        4.85787082     0.4699%
18            4.9759        5.0070631     -0.6224%
19            5.1755        5.1521962      0.4523%
20            5.2792        5.29358459    -0.2717%

Mathematica / Wolfram Language[edit]

[email protected][
Table[{n, #[[1]], #[[2]],
Row[{Round[10000 Abs[#[[1]] - #[[2]]]/#[[2]]]/100., "%"}]} &@
N[{Mean[Array[
[email protected][#, 1, UnsameQ[##] &, All] - 1 &[# /.
MapIndexed[#2[[1]] -> #1 &,
RandomInteger[{1, n}, n]] &] &, 10000]],
Sum[n! n^(n - k - 1)/(n - k)!, {k, n}]/n^(n - 1)}, 5], {n, 1,
20}], {"N", "average", "analytical", "error"}]
Output:
N	average	analytical	error
1	1.0000	1.0000		0.%
2	1.5017	1.5000		0.11%
3	1.8910	1.8889		0.11%
4	2.2334	2.2188		0.66%
5	2.5090	2.5104		0.06%
6	2.8092	2.7747		1.24%
7	3.0468	3.0181		0.95%
8	3.2253	3.2450		0.61%
9	3.4695	3.4583		0.32%
10	3.6661	3.6602		0.16%
11	3.8662	3.8524		0.36%
12	4.0393	4.0361		0.08%
13	4.2232	4.2123		0.26%
14	4.3496	4.3820		0.74%
15	4.5706	4.5458		0.55%
16	4.6963	4.7043		0.17%
17	4.8548	4.8579		0.06%
18	5.0671	5.0071		1.2%
19	5.1702	5.1522		0.35%
20	5.2264	5.2936		1.27%

Nim[edit]

Translation of: C
import math, strfmt
randomize()
 
const
maxN = 20
times = 1_000_000
 
proc factorial(n): float =
result = 1
for i in 1 .. n:
result *= i.float
 
proc expected(n): float =
for i in 1 .. n:
result += factorial(n) / pow(n.float, i.float) / factorial(n - i)
 
proc test(n, times): int =
for i in 1 .. times:
var
x = 1
bits = 0
while (bits and x) == 0:
inc result
bits = bits or x
x = 1 shl random(n)
 
echo " n\tavg\texp.\tdiff"
echo "-------------------------------"
for n in 1 .. maxN:
let cnt = test(n, times)
let avg = cnt.float / times
let theory = expected(n)
let diff = (avg / theory - 1) * 100
printlnfmt "{:2} {:8.4f} {:8.4f} {:6.3f}%", n, avg, theory, diff
Output:
 n	avg	exp.	diff
-------------------------------
 1   1.0000   1.0000      0%
 2   1.5001   1.5000  0.008%
 3   1.8884   1.8889 -0.025%
 4   2.2187   2.2187 -0.000%
 5   2.5098   2.5104 -0.025%
 6   2.7752   2.7747  0.017%
 7   3.0175   3.0181 -0.020%
 8   3.2411   3.2450 -0.120%
 9   3.4565   3.4583 -0.054%
10   3.6599   3.6602 -0.010%
11   3.8555   3.8524  0.081%
12   4.0381   4.0361  0.051%
13   4.2124   4.2123  0.000%
14   4.3813   4.3820 -0.017%
15   4.5471   4.5458  0.027%
16   4.7009   4.7043 -0.072%
17   4.8589   4.8579  0.021%
18   5.0054   5.0071 -0.034%
19   5.1554   5.1522  0.061%
20   5.2915   5.2936 -0.040%

PARI/GP[edit]

Translation of: C
expected(n)=sum(i=1,n,n!/(n-i)!/n^i,0.);
test(n, times)={
my(ct);
for(i=1,times,
my(x=1,bits);
while(!bitand(bits,x),ct++; bits=bitor(bits,x); x = 1<<random(n))
);
ct
};
TIMES=1000000;
{for(n=1,20,
my(cnt=test(n, TIMES),avg=cnt/TIMES,ex=expected(n),diff=(avg/ex-1)*100.);
print(n"\t"avg*1."\t"ex*1."\t"diff);
)}
Output:
1       1.0000  1.0000  0.E-7
2       1.4998  1.5000  -0.012933
3       1.8891  1.8889  0.013559
4       2.2198  2.2188  0.047369
5       2.5095  2.5104  -0.034616
6       2.7744  2.7747  -0.010248
7       3.0177  3.0181  -0.012945
8       3.2467  3.2450  0.050600
9       3.4611  3.4583  0.080278
10      3.6595  3.6602  -0.018651
11      3.8541  3.8524  0.044880
12      4.0428  4.0361  0.16690
13      4.2116  4.2123  -0.017921
14      4.3825  4.3820  0.011150
15      4.5467  4.5458  0.020562
16      4.7087  4.7043  0.095058
17      4.8573  4.8579  -0.011997
18      5.0080  5.0071  0.018312
19      5.1530  5.1522  0.015970
20      5.2970  5.2936  0.065143

Perl 6[edit]

Works with: rakudo version 2016.08
constant MAX_N  = 20;
constant TRIALS = 100;
 
for 1 .. MAX_N -> $N {
my $empiric = TRIALS R/ [+] find-loop(random-mapping($N)).elems xx TRIALS;
my $theoric = [+]
map -> $k { $N ** ($k + 1) R/ [*] flat $k**2, $N - $k + 1 .. $N }, 1 .. $N;
 
FIRST say " N empiric theoric (error)";
FIRST say "=== ========= ============ =========";
 
printf "%3d  %9.4f  %12.4f (%4.2f%%)\n",
$N, $empiric,
$theoric, 100 * abs($theoric - $empiric) / $theoric;
}
 
sub random-mapping { hash .list Z=> .roll given ^$^size }
sub find-loop { 0, | %^mapping{*} ...^ { (%){$_}++ } }
Example:
 N    empiric      theoric      (error)
===  =========  ============  =========
  1     1.0000        1.0000    (0.00%)
  2     1.5600        1.5000    (4.00%)
  3     1.7800        1.8889    (5.76%)
  4     2.1800        2.2188    (1.75%)
  5     2.6200        2.5104    (4.37%)
  6     2.8300        2.7747    (1.99%)
  7     3.1200        3.0181    (3.37%)
  8     3.1400        3.2450    (3.24%)
  9     3.4500        3.4583    (0.24%)
 10     3.6700        3.6602    (0.27%)
 11     3.8300        3.8524    (0.58%)
 12     4.3600        4.0361    (8.03%)
 13     3.9000        4.2123    (7.42%)
 14     4.4900        4.3820    (2.46%)
 15     4.9500        4.5458    (8.89%)
 16     4.9800        4.7043    (5.86%)
 17     4.9100        4.8579    (1.07%)
 18     4.9700        5.0071    (0.74%)
 19     5.1000        5.1522    (1.01%)
 20     5.2300        5.2936    (1.20%)

Phix[edit]

constant MAX = 20,
ITER = 1000000
 
function expected(integer n)
atom sum = 0
for i=1 to n do
sum += factorial(n) / power(n,i) / factorial(n-i)
end for
return sum
end function
 
function test(integer n)
integer count = 0, x, bits
for i=1 to ITER do
x = 1
bits = 0
while not and_bits(bits,x) do
count += 1
bits = or_bits(bits,x)
x = power(2,rand(n)-1)
end while
end for
return count/ITER
end function
 
atom av, ex
puts(1," n avg. exp. (error%)\n");
puts(1,"== ====== ====== ========\n");
for n=1 to MAX do
av = test(n)
ex = expected(n)
printf(1,"%2d %8.4f %8.4f (%5.3f%%)\n", {n,av,ex,abs(1-av/ex)*100})
end for
Output:
 n     avg.     exp.  (error%)
==   ======   ======  ========
 1   1.0000   1.0000  (0.000%)
 2   1.5003   1.5000  (0.018%)
 3   1.8880   1.8889  (0.046%)
 4   2.2176   2.2188  (0.052%)
 5   2.5104   2.5104  (0.001%)
 6   2.7734   2.7747  (0.046%)
 7   3.0198   3.0181  (0.055%)
 8   3.2464   3.2450  (0.042%)
 9   3.4562   3.4583  (0.062%)
10   3.6618   3.6602  (0.043%)
11   3.8511   3.8524  (0.033%)
12   4.0357   4.0361  (0.009%)
13   4.2158   4.2123  (0.083%)
14   4.3843   4.3820  (0.052%)
15   4.5410   4.5458  (0.105%)
16   4.7084   4.7043  (0.087%)
17   4.8603   4.8579  (0.049%)
18   5.0044   5.0071  (0.052%)
19   5.1516   5.1522  (0.011%)
20   5.2955   5.2936  (0.037%)

PicoLisp[edit]

Translation of: Python
(scl 4)
(seed (in "/dev/urandom" (rd 8)))
 
(de fact (N)
(if (=0 N) 1 (apply * (range 1 N))) )
 
(de analytical (N)
(sum
'((I)
(/
(* (fact N) 1.0)
(** N I)
(fact (- N I)) ) )
(range 1 N) ) )
 
(de testing (N)
(let (C 0 N (dec N) X 0 B 0 I 1000000)
(do I
(zero B)
(one X)
(while (=0 (& B X))
(inc 'C)
(setq
B (| B X)
X (** 2 (rand 0 N)) ) ) )
(*/ C 1.0 I) ) )
 
(let F (2 8 8 6)
(tab F "N" "Avg" "Exp" "Diff")
(for I 20
(let (A (testing I) B (analytical I))
(tab F
I
(round A 4)
(round B 4)
(round
(*
(abs (- (*/ A 1.0 B) 1.0))
100 )
2 ) ) ) ) )
 
(bye)

PowerShell[edit]

Works with: PowerShell version 2
 
function Get-AnalyticalLoopAverage ( [int]$N )
{
# Expected loop average = sum from i = 1 to N of N! / (N-i)! / N^(N-i+1)
# Equivalently, Expected loop average = sum from i = 1 to N of F(i)
# where F(N) = 1, and F(i) = F(i+1)*i/N
 
$LoopAverage = $Fi = 1
 
If ( $N -eq 1 ) { return $LoopAverage }
 
ForEach ( $i in ($N-1)..1 )
{
$Fi *= $i / $N
$LoopAverage += $Fi
}
return $LoopAverage
}
 
function Get-ExperimentalLoopAverage ( [int]$N, [int]$Tests = 100000 )
{
If ( $N -eq 1 ) { return 1 }
 
# Using 0 through N-1 instead of 1 through N for speed and simplicity
$NMO = $N - 1
 
# Create array to hold mapping function
$F = New-Object int[] ( $N )
 
$Count = 0
$Random = New-Object System.Random
 
ForEach ( $Test in 1..$Tests )
{
# Map each number to a random number
ForEach ( $i in 0..$NMO )
{
$F[$i] = $Random.Next( $N )
}
 
# For each number...
ForEach ( $i in 0..$NMO )
{
# Add the number to the list
$List = @()
$Count++
$List += $X = $i
 
# If loop does not yet exist in list...
While ( $F[$X] -notin $List )
{
# Go to the next mapped number and add it to the list
$Count++
$List += $X = $F[$X]
}
}
}
$LoopAvereage = $Count / $N / $Tests
return $LoopAvereage
}
 

Note: The use of the [pscustomobject] type accelerator to simplify making the test result table look pretty requires PowerShell 3.0.

 
# Display results for N = 1 through 20
ForEach ( $N in 1..20 )
{
$AnalyticalAverage = Get-AnalyticalLoopAverage $N
$ExperimentalAverage = Get-ExperimentalLoopAverage $N
[pscustomobject] @{
N = $N.ToString().PadLeft( 2, ' ' )
Analytical = $AnalyticalAverage.ToString( '0.00000000' )
Experimental = $ExperimentalAverage.ToString( '0.00000000' )
'Error (%)' = ( [math]::Abs( $AnalyticalAverage - $ExperimentalAverage ) / $AnalyticalAverage * 100 ).ToString( '0.00000000' )
}
}
 
Output:
N  Analytical Experimental Error (%) 
-  ---------- ------------ --------- 
 1 1.00000000 1.00000000   0.00000000
 2 1.50000000 1.49985500   0.00966667
 3 1.88888889 1.88713000   0.09311765
 4 2.21875000 2.22103500   0.10298592
 5 2.51040000 2.51069200   0.01163161
 6 2.77469136 2.77264833   0.07363070
 7 3.01813870 3.01547143   0.08837474
 8 3.24501801 3.25003875   0.15472163
 9 3.45831574 3.45067667   0.22089013
10 3.66021568 3.65659000   0.09905646
11 3.85237205 3.85669273   0.11215626
12 4.03607368 4.03813500   0.05107253
13 4.21234791 4.20946231   0.06850349
14 4.38202942 4.38458786   0.05838465
15 4.54580729 4.54466400   0.02515032
16 4.70425825 4.70146375   0.05940356
17 4.85787082 4.86807647   0.21008483
18 5.00706310 5.01939278   0.24624572
19 5.15219620 5.15179263   0.00783296
20 5.29358459 5.29214950   0.02710991

Python[edit]

Translation of: C
from __future__ import division # Only necessary for Python 2.X
from math import factorial
from random import randrange
 
MAX_N = 20
TIMES = 1000000
 
def analytical(n):
return sum(factorial(n) / pow(n, i) / factorial(n -i) for i in range(1, n+1))
 
def test(n, times):
count = 0
for i in range(times):
x, bits = 1, 0
while not (bits & x):
count += 1
bits |= x
x = 1 << randrange(n)
return count / times
 
if __name__ == '__main__':
print(" n\tavg\texp.\tdiff\n-------------------------------")
for n in range(1, MAX_N+1):
avg = test(n, TIMES)
theory = analytical(n)
diff = (avg / theory - 1) * 100
print("%2d %8.4f %8.4f %6.3f%%" % (n, avg, theory, diff))
Output:
 n	avg	exp.	diff
-------------------------------
 1   1.0000   1.0000  0.000%
 2   1.5006   1.5000  0.037%
 3   1.8887   1.8889 -0.012%
 4   2.2190   2.2188  0.011%
 5   2.5101   2.5104 -0.012%
 6   2.7750   2.7747  0.012%
 7   3.0158   3.0181 -0.076%
 8   3.2447   3.2450 -0.009%
 9   3.4586   3.4583  0.009%
10   3.6598   3.6602 -0.010%
11   3.8510   3.8524 -0.036%
12   4.0368   4.0361  0.017%
13   4.2099   4.2123 -0.058%
14   4.3784   4.3820 -0.083%
15   4.5484   4.5458  0.058%
16   4.7045   4.7043  0.006%
17   4.8611   4.8579  0.067%
18   5.0074   5.0071  0.007%
19   5.1534   5.1522  0.024%
20   5.2927   5.2936 -0.017%

Racket[edit]

 
#lang racket
(require (only-in math factorial))
 
(define (analytical n)
(for/sum ([i (in-range 1 (add1 n))])
(/ (factorial n) (expt n i) (factorial (- n i)))))
 
(define (test n times)
(define (count-times seen times)
(define x (random n))
(if (memq x seen) times (count-times (cons x seen) (add1 times))))
(/ (for/fold ([count 0]) ([i times]) (count-times '() count))
times))
 
(define (test-table max-n times)
(displayln " n avg theory error\n------------------------")
(for ([i (in-range 1 (add1 max-n))])
(define average (test i times))
(define theory (analytical i))
(define difference (* (abs (sub1 (/ average theory))) 100))
(displayln (~a (~a i #:width 2 #:align 'right)
" " (real->decimal-string average 4)
" " (real->decimal-string theory 4)
" " (real->decimal-string difference 4)
"%"))))
 
(test-table 20 10000)
 
Output:
 n avg    theory error
------------------------
 1 1.0000 1.0000 0.0000%
 2 1.5082 1.5000 0.5467%
 3 1.8966 1.8889 0.4082%
 4 2.2251 2.2188 0.2862%
 5 2.5138 2.5104 0.1354%
 6 2.7582 2.7747 0.5943%
 7 3.0253 3.0181 0.2373%
 8 3.2293 3.2450 0.4844%
 9 3.4602 3.4583 0.0545%
10 3.6831 3.6602 0.6252%
11 3.8459 3.8524 0.1680%
12 4.0348 4.0361 0.0316%
13 4.1896 4.2123 0.5400%
14 4.3555 4.3820 0.6054%
15 4.5678 4.5458 0.4838%
16 4.6950 4.7043 0.1968%
17 4.8524 4.8579 0.1126%
18 5.0224 5.0071 0.3063%
19 5.1017 5.1522 0.9801%
20 5.3316 5.2936 0.7181%

REXX[edit]

This REXX program automatically adjusts the precision (decimal digits) to be used based on the size of the
factorial (product) for   RUNS.

Also note that the   !   (factorial function)   uses memoization for optimization.

/*REXX program computes the average loop length mapping a random field 1···N ───► 1···N */
parse arg runs tests seed . /*obtain optional arguments from the CL*/
if runs =='' | runs =="," then runs = 40 /*Not specified? Then use the default.*/
if tests =='' | tests =="," then tests= 1000000 /* " " " " " " */
if datatype(seed,'W') then call random ,, seed /*Is integer? For RAND repeatability.*/
!.=0;  !.0=1 /*used for factorial (!) memoization.*/
numeric digits 100000 /*be able to calculate 25k! if need be.*/
numeric digits max(9, length( !(runs) ) ) /*set the NUMERIC DIGITS for  !(runs). */
say right( runs, 24) 'runs' /*display number of runs we're using.*/
say right( tests, 24) 'tests' /* " " " tests " " */
say right( digits(), 24) 'digits' /* " " " digits " " */
say
say " N average exact  % error " /* ◄─── title, header ►────────┐ */
hdr=" ═══ ═════════ ═════════ ═════════"; pad=left('',3) /* ◄────────┘ */
say hdr
do #=1 for runs; av=fmtD(exact(#)) /*use four digits past decimal point. */
xa=fmtD(exper(#)) /* " " " " " " */
say right(#,9) pad xa pad av pad fmtD(abs(xa-av)*100/av) /*display values.*/
end /*#*/
say hdr /*display the final header (some bars).*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
!: procedure expose !.; parse arg z; if !.z\==0 then return !.z
 !=1; do j=2 for z-1;  !=!*j;  !.j=!; end; /*factorial*/ return !
/*──────────────────────────────────────────────────────────────────────────────────────*/
exact: parse arg x; s=0; do j=1 for x; s=s+!(x)/!(x-j)/x**j; end; return s
/*──────────────────────────────────────────────────────────────────────────────────────*/
exper: parse arg n; k=0; do tests; $.=0 /*do it TESTS times.*/
do n; r=random(1,n); if $.r then leave
$.r=1; k=k+1 /*bump the counter. */
end /*n*/
end /*tests*/
return k/tests
/*──────────────────────────────────────────────────────────────────────────────────────*/
fmtD: parse arg y,d; d=word(d 4,1); y=format(y,,d); parse var y w '.' f
if f=0 then return w || left('', d+1); return y

output   when using the default inputs:

                      40 runs
                 1000000 tests
                      48 digits

         N    average     exact     % error
        ═══  ═════════  ═════════  ═════════
         1     1          1          0
         2     1.4964     1.5000     0.2400
         3     1.8876     1.8889     0.0688
         4     2.2222     2.2188     0.1532
         5     2.5104     2.5104     0
         6     2.7758     2.7747     0.0396
         7     3.0194     3.0181     0.0431
         8     3.2608     3.2450     0.4869
         9     3.4565     3.4583     0.0520
        10     3.6583     3.6602     0.0519
        11     3.8513     3.8524     0.0286
        12     4.0401     4.0361     0.0991
        13     4.2133     4.2123     0.0237
        14     4.3835     4.3820     0.0342
        15     4.5445     4.5458     0.0286
        16     4.6672     4.7043     0.7886
        17     4.8575     4.8579     0.0082
        18     5.0105     5.0071     0.0679
        19     5.1517     5.1522     0.0097
        20     5.2903     5.2936     0.0623
        21     5.4328     5.4315     0.0239
        22     5.5674     5.5662     0.0216
        23     5.6990     5.6979     0.0193
        24     5.8353     5.8268     0.1459
        25     5.9536     5.9530     0.0101
        26     6.0801     6.0767     0.0560
        27     6.1997     6.1981     0.0258
        28     6.3197     6.3173     0.0380
        29     6.4328     6.4344     0.0249
        30     6.5485     6.5495     0.0153
        31     6.6615     6.6627     0.0180
        32     6.7102     6.7740     0.9418
        33     6.8826     6.8837     0.0160
        34     6.9878     6.9917     0.0558
        35     7.0996     7.0982     0.0197
        36     7.2054     7.2032     0.0305
        37     7.3073     7.3067     0.0082
        38     7.4089     7.4088     0.0013
        39     7.5052     7.5096     0.0586
        40     7.6151     7.6091     0.0789
        ═══  ═════════  ═════════  ═════════ 

Ruby[edit]

Ruby does not have a factorial method, not even in it's math library.

class Integer
def factorial
self == 0 ? 1 : (1..self).inject(:*)
end
end
 
def rand_until_rep(n)
rands = {}
loop do
r = rand(1..n)
return rands.size if rands[r]
rands[r] = true
end
end
 
runs = 1_000_000
 
puts " N average exp. diff ",
"=== ======== ======== ==========="
(1..20).each do |n|
sum_of_runs = runs.times.inject(0){|sum, _| sum += rand_until_rep(n)}
avg = sum_of_runs / runs.to_f
analytical = (1..n).inject(0){|sum, i| sum += (n.factorial / (n**i).to_f / (n-i).factorial)}
puts "%3d  %8.4f  %8.4f (%8.4f%%)" % [n, avg, analytical, (avg/analytical - 1)*100]
end
Output:
 N    average    exp.        diff   
===  ========  ========  ===========
  1    1.0000    1.0000  (  0.0000%)
  2    1.4999    1.5000  ( -0.0054%)
  3    1.8886    1.8889  ( -0.0158%)
  4    2.2181    2.2188  ( -0.0293%)
  5    2.5107    2.5104  (  0.0110%)
  6    2.7717    2.7747  ( -0.1074%)
  7    3.0167    3.0181  ( -0.0484%)
  8    3.2442    3.2450  ( -0.0257%)
  9    3.4597    3.4583  (  0.0394%)
 10    3.6572    3.6602  ( -0.0821%)
 11    3.8502    3.8524  ( -0.0562%)
 12    4.0357    4.0361  ( -0.0084%)
 13    4.2139    4.2123  (  0.0360%)
 14    4.3805    4.3820  ( -0.0360%)
 15    4.5481    4.5458  (  0.0505%)
 16    4.7030    4.7043  ( -0.0265%)
 17    4.8582    4.8579  (  0.0075%)
 18    5.0078    5.0071  (  0.0151%)
 19    5.1568    5.1522  (  0.0893%)
 20    5.2885    5.2936  ( -0.0961%)

Rust[edit]

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
}
 
 
 
Output:

Using default arguments:

N    average    analytical    (error)
===  =========  ============  =========
  1     1.0000        1.0000  ( +0.00%)
  2     1.4992        1.5000  ( -0.05%)
  3     1.8881        1.8889  ( -0.04%)
  4     2.2177        2.2188  ( -0.05%)
  5     2.5107        2.5104  ( +0.01%)
  6     2.7752        2.7747  ( +0.02%)
  7     3.0172        3.0181  ( -0.03%)
  8     3.2452        3.2450  ( +0.01%)
  9     3.4628        3.4583  ( +0.13%)
 10     3.6606        3.6602  ( +0.01%)
 11     3.8515        3.8524  ( -0.02%)
 12     4.0348        4.0361  ( -0.03%)
 13     4.2105        4.2123  ( -0.04%)
 14     4.3835        4.3820  ( +0.03%)
 15     4.5477        4.5458  ( +0.04%)
 16     4.7042        4.7043  ( -0.00%)
 17     4.8580        4.8579  ( +0.00%)
 18     5.0076        5.0071  ( +0.01%)
 19     5.1554        5.1522  ( +0.06%)
 20     5.2911        5.2936  ( -0.05%)

Scala[edit]

 
import scala.util.Random
 
object AverageLoopLength extends App {
 
val factorial: Stream[Double] = 1 #:: factorial.zip(Stream.from(1)).map(n => n._2 * factorial(n._2 - 1))
 
def expected(n: Int) = (for (i <- 1 to n) yield factorial(n) / Math.pow(n, i) / factorial(n - i)).sum
 
def trial(n: Int):Double = {
var count = 0
var x = 1
var bits = 0
 
while ((bits & x) == 0) {
count = count + 1
bits = bits | x
x = 1 << Random.nextInt(n)
}
count
}
 
def tested(n: Int, times: Int) = (for (i <- 1 to times) yield trial(n)).sum / times
 
val results = for (n <- 1 to 20;
avg = tested(n, 1000000);
theory = expected(n)
) yield (n, avg, theory, (avg / theory - 1) * 100)
 
 
println("n avg exp diff")
println("------------------------------------")
results foreach { n => {
println(f"${n._1}%2d ${n._2}%2.6f ${n._3}%2.6f ${n._4}%2.3f%%")
}
}
 
}
 
Output:
n          avg         exp      diff
------------------------------------
 1    1.000000    1.000000    0.000%
 2    1.499894    1.500000    -0.007%
 3    1.887826    1.888889    -0.056%
 4    2.217514    2.218750    -0.056%
 5    2.510049    2.510400    -0.014%
 6    2.773658    2.774691    -0.037%
 7    3.016585    3.018139    -0.051%
 8    3.246865    3.245018    0.057%
 9    3.458683    3.458316    0.011%
10    3.660361    3.660216    0.004%
11    3.852663    3.852372    0.008%
12    4.036970    4.036074    0.022%
13    4.213653    4.212348    0.031%
14    4.385226    4.382029    0.073%
15    4.545667    4.545807    -0.003%
16    4.705559    4.704258    0.028%
17    4.854056    4.857871    -0.079%
18    5.007146    5.007063    0.002%
19    5.148767    5.152196    -0.067%
20    5.292875    5.293585    -0.013%

Scheme[edit]

 
(import (scheme base)
(scheme write)
(srfi 1 lists)
(only (srfi 13 strings) string-pad-right)
(srfi 27 random-bits))
 
(define (analytical-function n)
(define (factorial n)
(fold * 1 (iota n 1)))
;
(fold (lambda (i sum)
(+ sum
(/ (factorial n) (expt n i) (factorial (- n i)))))
0
(iota n 1)))
 
(define (simulation n runs)
(define (single-simulation)
(random-source-randomize! default-random-source)
(let ((vec (make-vector n #f)))
(let loop ((count 0)
(num (random-integer n)))
(if (vector-ref vec num)
count
(begin (vector-set! vec num #t)
(loop (+ 1 count)
(random-integer n)))))))
;;
(let loop ((total 0)
(run runs))
(if (zero? run)
(/ total runs)
(loop (+ total (single-simulation))
(- run 1)))))
 
(display " N average formula (error) \n")
(display "=== ========= ========= =========\n")
(for-each
(lambda (n)
(let ((simulation (inexact (simulation n 10000)))
(formula (inexact (analytical-function n))))
(display
(string-append
" "
(string-pad-right (number->string n) 3)
" "
(string-pad-right (number->string simulation) 6)
" "
(string-pad-right (number->string formula) 6)
" ("
(string-pad-right
(number->string (* 100 (/ (- simulation formula) formula)))
5)
"%)"))
(newline)))
(iota 20 1))
 
Output:
 N   average   formula   (error) 
=== ========= ========= =========
 1     1.0      1.0      (0.0  %)
 2     1.5018   1.5      (0.120%)
 3     1.8863   1.8888   (-0.13%)
 4     2.2154   2.2187   (-0.15%)
 5     2.5082   2.5104   (-0.08%)
 6     2.7613   2.7746   (-0.48%)
 7     3.036    3.0181   (0.591%)
 8     3.2656   3.2450   (0.634%)
 9     3.455    3.4583   (-0.09%)
 10    3.682    3.6602   (0.595%)
 11    3.8233   3.8523   (-0.75%)
 12    4.0409   4.0360   (0.119%)
 13    4.2471   4.2123   (0.825%)
 14    4.3577   4.3820   (-0.55%)
 15    4.5351   4.5458   (-0.23%)
 16    4.7181   4.7042   (0.294%)
 17    4.8877   4.8578   (0.614%)
 18    5.0239   5.0070   (0.336%)
 19    5.1216   5.1521   (-0.59%)
 20    5.2717   5.2935   (-0.41%)

Seed7[edit]

$ include "seed7_05.s7i";
include "float.s7i";
 
const integer: TESTS is 1000000;
 
const func float: factorial (in integer: number) is func
result
var float: factorial is 1.0;
local
var integer: i is 0;
begin
for i range 2 to number do
factorial *:= flt(i);
end for;
end func;
 
const func float: analytical (in integer: number) is func
result
var float: sum is 0.0;
local
var integer: i is 0;
begin
for i range 1 to number do
sum +:= factorial(number) / factorial(number - i) / flt(number)**i;
end for;
end func;
 
const func float: experimental (in integer: number) is func
result
var float: experimental is 0.0;
local
var integer: run is 0;
var set of integer: seen is EMPTY_SET;
var integer: current is 1;
var integer: count is 0;
begin
for run range 1 to TESTS do
current := 1;
seen := EMPTY_SET;
while current not in seen do
incr(count);
incl(seen, current);
current := rand(1, number);
end while;
end for;
experimental := flt(count) / flt(TESTS);
end func;
 
const proc: main is func
local
var integer: number is 0;
var float: analytical is 0.0;
var float: experimental is 0.0;
var float: err is 0.0;
begin
writeln(" N avg calc  %diff");
for number range 1 to 20 do
analytical := analytical(number);
experimental := experimental(number);
err := abs(experimental - analytical) / analytical * 100.0;
writeln(number lpad 2 <& experimental digits 4 lpad 7 <&
analytical digits 4 lpad 7 <& err digits 3 lpad 7);
end for;
end func;
Output:
 N  avg    calc   %diff
 1 1.0000 1.0000  0.000
 2 1.4999 1.5000  0.005
 3 1.8891 1.8889  0.010
 4 2.2196 2.2188  0.040
 5 2.5073 2.5104  0.122
 6 2.7744 2.7747  0.010
 7 3.0186 3.0181  0.015
 8 3.2463 3.2450  0.040
 9 3.4592 3.4583  0.027
10 3.6597 3.6602  0.013
11 3.8549 3.8524  0.066
12 4.0374 4.0361  0.033
13 4.2115 4.2123  0.019
14 4.3835 4.3820  0.033
15 4.5474 4.5458  0.035
16 4.7017 4.7043  0.055
17 4.8558 4.8579  0.043
18 5.0096 5.0071  0.051
19 5.1522 5.1522  0.000
20 5.2907 5.2936  0.054

Tcl[edit]

# Generate a list of the numbers increasing from $a to $b
proc range {a b} {
for {set result {}} {$a <= $b} {incr a} {lappend result $a}
return $result
}
 
# Computing the expected value analytically
proc tcl::mathfunc::factorial n {
 ::tcl::mathop::* {*}[range 2 $n]
}
proc Analytical {n} {
set sum 0.0
foreach x [range 1 $n] {
set sum [expr {$sum + factorial($n) / factorial($n-$x) / double($n)**$x}]
}
return $sum
}
 
# Determining an approximation to the value experimentally
proc Experimental {n numTests} {
set count 0
set u0 [lrepeat $n 1]
foreach run [range 1 $numTests] {
set unseen $u0
for {set i 0} {[lindex $unseen $i]} {incr count} {
lset unseen $i 0
set i [expr {int(rand()*$n)}]
}
}
return [expr {$count / double($numTests)}]
}
 
# Tabulate the results in exactly the original format
puts " N average analytical (error)"
puts "=== ========= ============ ========="
foreach n [range 1 20] {
set a [Analytical $n]
set e [Experimental $n 100000]
puts [format "%3d  %9.4f  %12.4f (%6.2f%%)" $n $e $a [expr {abs($e-$a)/$a*100.0}]]
}
Output:
 N    average    analytical    (error)
===  =========  ============  =========
  1     1.0000        1.0000  (  0.00%)
  2     1.5003        1.5000  (  0.02%)
  3     1.8881        1.8889  (  0.04%)
  4     2.2228        2.2188  (  0.18%)
  5     2.5109        2.5104  (  0.02%)
  6     2.7804        2.7747  (  0.20%)
  7     3.0223        3.0181  (  0.14%)
  8     3.2456        3.2450  (  0.02%)
  9     3.4598        3.4583  (  0.04%)
 10     3.6590        3.6602  (  0.03%)
 11     3.8527        3.8524  (  0.01%)
 12     4.0390        4.0361  (  0.07%)
 13     4.2156        4.2123  (  0.08%)
 14     4.3821        4.3820  (  0.00%)
 15     4.5527        4.5458  (  0.15%)
 16     4.6952        4.7043  (  0.19%)
 17     4.8530        4.8579  (  0.10%)
 18     4.9912        5.0071  (  0.32%)
 19     5.1578        5.1522  (  0.11%)
 20     5.2992        5.2936  (  0.11%)

Unicon[edit]

Translation of: C
link printf, factors
 
$define MAX_N 20
$define TIMES 1000000
$define RAND_MAX 2147483647
 
procedure expected(n)
local sum := 0
every i := 1 to n do
sum +:= factorial(n) / (n ^ i) / factorial(n - i)
return sum
end
 
procedure test(n, times)
local i, count := 0, x, bits
every i := 0 to times-1 do {
x := 1
bits := 0
while iand(bits, x)=0 do {
count +:= 1
bits := ior(bits, x)
x := ishift(1 , ?n-1)
}
}
return count
end
 
procedure main(void)
local n, cnt, avg, theory, diff
write(" n\tavg\texp.\tdiff\n", repl("-",29))
every n := 1 to MAX_N do {
cnt := test(n, TIMES)
avg := real(cnt) / TIMES
theory := expected(n)
diff := (avg / theory - 1) * 100
printf("%2d %8.4r %8.4r %6.3r%%\n", n, avg, theory, diff)
}
return 0
end
Output:
 n      avg     exp.    diff                                                    
-----------------------------                                                   
 1   1.0000   1.0000  0.000%                                                    
 2   1.5008   1.5000  0.056%                                                    
 3   1.8879   1.8889 -0.051%                                                    
 4   2.2208   2.2188  0.091%                                                    
 5   2.5127   2.5104  0.093%                                                    
 6   2.7759   2.7747  0.044%                                                    
 7   3.0175   3.0181 -0.023%                                                    
 8   3.2425   3.2450 -0.079%                                                    
 9   3.4571   3.4583 -0.034%                                                    
10   3.6613   3.6602  0.029%                                                    
11   3.8493   3.8524 -0.081%                                                    
12   4.0384   4.0361  0.058%                                                    
13   4.2133   4.2123  0.023%                                                    
14   4.3804   4.3820 -0.037%                                                    
15   4.5475   4.5458  0.038%                                                    
16   4.7049   4.7043  0.014%                                                    
17   4.8575   4.8579 -0.008%
18   5.0088   5.0071  0.035%
19   5.1533   5.1522  0.021%
20   5.2893   5.2936 -0.081%

zkl[edit]

const N=20;
 
(" N average analytical (error)").println();
("=== ========= ============ =========").println();
foreach n in ([1..N]){
a := avg(n);
b := ana(n);
"%3d  %9.4f  %12.4f (%6.2f%%)".fmt(
n, a, b, ((a-b)/b*100)).println();
}
 
fcn f(n){ (0).random(n) }
 
fcn avg(n){
tests := 0d10_000;
sum := 0;
do(tests){
v:=(0).pump(n,List,T(Void,False)).copy();
while(1){
z := f(n);
if(v[z]) break;
v[z] = True;
sum += 1;
}
}
return(sum.toFloat() / tests);
}
 
fcn fact(n) { (1).reduce(n,fcn(N,n){N*n},1.0) } //-->Float
fcn ana(n){
n=n.toFloat();
(1).reduce(n,'wrap(sum,i){ sum+fact(n)/n.pow(i)/fact(n-i) },0.0);
}
Output:
 N    average    analytical    (error)
===  =========  ============  =========
  1     1.0000        1.0000  (  0.00%)
  2     1.5053        1.5000  (  0.35%)
  3     1.8899        1.8889  (  0.05%)
  4     2.2384        2.2188  (  0.89%)
  5     2.5090        2.5104  ( -0.06%)
  6     2.7824        2.7747  (  0.28%)
  7     3.0449        3.0181  (  0.89%)
  8     3.2430        3.2450  ( -0.06%)
  9     3.4744        3.4583  (  0.47%)
 10     3.6693        3.6602  (  0.25%)
 11     3.8833        3.8524  (  0.80%)
 12     4.0225        4.0361  ( -0.34%)
 13     4.1899        4.2123  ( -0.53%)
 14     4.4135        4.3820  (  0.72%)
 15     4.5807        4.5458  (  0.77%)
 16     4.7304        4.7043  (  0.56%)
 17     4.8437        4.8579  ( -0.29%)
 18     4.9838        5.0071  ( -0.46%)
 19     5.1767        5.1522  (  0.48%)
 20     5.2723        5.2936  ( -0.40%)