Given a probability $p$ and an integer $N$, I would like to generate a sample $S$ of the population $P=\{0,1,...,N\}$ such that integer $m\in P$ is sampled with probability $p^m$.
It is trivial to do this in time $O(N)$. I am looking for an algorithm that can perform this sampling in (expected) time proportional to the output, say $O(|S|)=O(\tfrac{1-p^{N+1}}{1-p})$.
Numerical precision is very important, and in most cases $p$ is very close to $1$ so that the expected size of the sample will remain constant as $N$ increases.
If the sampling probability were constant I could use a Poisson distribution to "jump" over unsampled integers (among other methods), but in this case the probability is not constant. What do I do?
