function n = geo2rnd(a, b)
r = rand;
n = 1;
prb = a;
while r > prb
    n = n + 1;
    prb = prb + (1-b)^(n-2)*(1-a)*b;
end
n = 2*n;