F ffactorial(n)
   V result = 1.0
   L(i) 2..n
      result *= i
   R result

V MAX_N = 20
V TIMES = 1000000

F analytical(n)
   R sum((1..n).map(i -> ffactorial(@n) / pow(Float(@n), Float(i)) / ffactorial(@n - i)))

F test(n, times)
   V count = 0
   L(i) 0 .< times
      V (x, bits) = (1, 0)
      L (bits [&] x) == 0
         count++
         bits [|]= x
         x = 1 << random:(n)
   R Float(count) / times

print(" n      avg     exp.    diff\n-------------------------------")
L(n) 1 .. MAX_N
   V avg = test(n, TIMES)
   V theory = analytical(n)
   V diff = (avg / theory - 1) * 100
   print(‘#2 #3.4 #3.4 #2.3%’.format(n, avg, theory, diff))