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