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