include xpllib;         \for Print
def MAX_N = 20;
def TIMES = 1000000;

func real Factorial_(N);
int  N;
real F;
int  I;
[F:= 1.;
for I:= 1 to N do F:= F * float(I);
return F;
];

func real Expected(N);
int  N;
real Sum;
int  I;
[Sum:= 0.;
for I:= 1 to N do
    Sum:= Sum + Factorial_(N) / Pow(float(N), float(I)) / Factorial_(N-I);
return Sum;
];

func RandInt(N);
int  N;
int  R, RMax;
def  RAND_MAX = -1>>1;
[RMax:= RAND_MAX / N * N;
repeat  R:= Ran(RAND_MAX);
until   R < RMax;
return R / (RAND_MAX / N);
];

func Test(N, Times);
int  N, Times;
int  I, Count, X, Bits;
[Count:= 0;
for I:= 0 to Times-1 do
    [X:= 1;  Bits:= 0;
    while (Bits & X) = 0 do
        [Count:= Count+1;
        Bits:= Bits ! X;
        X:= 1 << RandInt(N);
        ];
    ];
return Count;
];

int  N, Cnt;
real Avg, Theory, Diff;
[\\srand(time(0));      random seeding is automatically done
Print(" N\tAvg\tExp.\tDiff\n-------------------------------\n");
for N:= 1 to MAX_N do
    [Cnt:= Test(N, TIMES);
    Avg:= float(Cnt) / float(TIMES);
    Theory:= Expected(N);
    Diff:= (Avg / Theory - 1.) * 100.;
    Print("%2d %3.4f %3.4f %2.3f%%\n", N, Avg, Theory, Diff);
    ];
]