function n = geosample(p)
r = rand;
n = 1;
prb = p^2;
while r > prb
    n = n + 1;
    prb = prb + n*(1-p)^(n-1)*p^2;
end
n = 2*n;