<p>Devroye's <em><a href="http://luc.devroye.org/rnbookindex.html" rel="nofollow noreferrer">Non-Uniform Random Variate Generation</a></em>, pp. 505 and 86, mentions an inversion by sequential search algorithm.</p>
<p>Based on that algorithm, if you know the <code>mean</code> is considerably less than 1, then if you generate a uniform random variate <code>u</code> in [0, 1], the Poisson variable will be 0 if <code>u &lt;= exp(-mean)</code>, and greater than 0 otherwise.</p>
<p>If the mean is low and you can tolerate an approximate distribution, then you can use the following approach (see Appendix A of &quot;<a href="https://arxiv.org/pdf/2004.00010.pdf" rel="nofollow noreferrer">The Discrete Gaussian for Differential Privacy</a>&quot;):</p>
<ol>
<li>Express <code>mean</code> in the form of a rational number, in the form <code>numer</code>/<code>denom</code>. For example, if <code>mean</code> is a fixed value, then <code>numer</code> and <code>denom</code> can be precalculated accordingly, such as at compile time.</li>
<li>Randomly generate a Bernoulli(<code>numer / denom</code>) number (generate 1 with probability <code>numer / denom</code> or 0 otherwise). If 1 was generated this way, repeat this step with Bernoulli(<code>numer / (denom * 2)</code>), Bernoulli(<code>numer / (denom * 3)</code>), and so on until 0 is generated this way. Generate these numbers using an algorithm that minimizes waste of bits, such as the one mentioned in Appendix B of Lumbroso's Fast Dice Roller paper (2013) or the &quot;ZeroToOne&quot; method modified from there and given in my section on <a href="https://github.com/peteroupc/peteroupc.github.io/blob/master/randomfunc.md#boolean-truefalse-conditions" rel="nofollow noreferrer">Boolean conditions</a>. See also <a href="https://stackoverflow.com/questions/60777414/uniformly-distributed-bit-sequence">this question</a>.</li>
<li>If step 2 produced an even number of ones, the Poisson variable is exactly 0.</li>
<li>If step 2 produced an odd number of ones, the Poisson variable is greater than 0, and a &quot;slower&quot; algorithm is necessary that samples only Poisson variables greater than 0.</li>
</ol>
<p>For example, say the mean is 1e-6 (1/1000000), Generate a Bernoulli(1/1000000) number, then Bernoulli(1/2000000), etc. Until you generate 0 this way. If an even number of ones were generated, then the Poisson variable is exactly 0. Otherwise, the Poisson variable is 1 or greater and a &quot;slower&quot; algorithm is necessary.</p>
<p>One example is the algorithm below, which is based on the one from pages 505 and 86, but only samples Poisson variables 1 or greater:</p>
<pre><code>METHOD Poisson1OrGreater(mean)
 sum=Math.exp(-mean)
 prod=sum
 u=RNDRANGE(sum, 1)
 i=0
 while i==0 or u&gt;sum
   prod*=mean/(i+1)
   sum+=prod
   i=i+1
 end
 return i
END METHOD
</code></pre>
<p>This method, though, is not very robust, especially since it uses numbers close to 1 (where the floating-point space is more sparse) rather than numbers close to 0.</p>
<hr />
<p>Note that the sum of <code>n</code> independent Poisson(<code>mean</code>) random variates is Poisson(<code>mean*n</code>) distributed (p. 501). Thus, the discussion above in this answer applies to a sum of <code>n</code> Poisson random variates as long as <code>n</code> times their mean remains small. For example, to generate a sum of 1000 Poisson random variates with a mean of 1e-6, simply generate a single Poisson random variate with a mean of 0.001. This will save considerably on calls to the pseudorandom number generator.</p>
<hr />
<p>There is another way to generate Poisson variates with low mean (1 or less). It is described by Duchon and Duvignau in &quot;Preserving the number of cycles of length k in a growing uniform permutation&quot;, Electronic Journal of Combinatorics 23(4), 2016.</p>
<p>First, generate a Poisson(1) random variate <code>x = Poisson1()</code> using the algorithm given below which uses only integer arithmetic (where <code>RNDINT(a)</code> generates a uniform random integer in [0, <code>a</code>]):</p>
<pre><code>METHOD Poisson1()
  ret=1; a=1; b=0
  while true // until this method returns
    j=RNDINT(a)
    if j&lt;a and j&lt;b: return ret
    if j==a: ret=ret+1
    else
      ret=ret-1; b=a+1
    end
    a=a+1
  end
END METHOD
</code></pre>
<p>Now let <code>mean</code> be the desired mean. Flip a coin <code>x</code> times, where the coin shows heads with probability equal to <code>mean</code>. (In other words, generate a binomial(<code>x</code>, <code>mean</code>) random variate.) The number of heads is a Poisson(<code>mean</code>) random variate.</p>
Performance for drawing numbers from Poisson distribution with low mean
<p>In order to draw random number from a Poisson distribution in C++, it is generally advised to use</p>

<pre><code>RNG_type rng;
std::poisson_distribution&lt;size_t&gt; d(1e-6);
auto r = d(rng);
</code></pre>

<p>At each call of the <code>std::poisson_distribution</code> object, an entire sequence of random bits is consumed (e.g. 32 bits with <a href="http://www.cplusplus.com/reference/random/mt19937/" rel="nofollow noreferrer">std::mt19937</a>, 64 bits for <a href="http://www.cplusplus.com/reference/random/mt19937_64/" rel="nofollow noreferrer">std::mt19937_64</a>). It strikes me that with such low mean (<code>mean = 1e-6</code>), the vast majority of times, only a few bits are enough to determine that the value to return is 0. The other bits could then be cached for later use.</p>

<p>Assuming that a sequence of bits set to true is associated to a high returned value from the Poisson distribution, when using a mean of <code>1e-6</code>, any sequence not starting with 19 trues necessarily returns a zero! Indeed, </p>

<pre><code>1 - 1/2^19 &lt; P(0, 1e-6) &lt; 1 - 1/2^20
</code></pre>

<p>, where <code>P(n, r)</code> denotes the probability of drawing <code>n</code> from a Poisson distribution with mean <code>r</code>. An algorithm that does not waste bits would use one bit half of the time, two bits a quarter of the times, three bits an eighth of the times, ....</p>

<p><strong>Is there an algorithm out there that can improve performance by consuming as few bits as possible when drawing Poisson numbers? Is there another way to improve performance compared to <code>std::poisson_distribution</code> when we consider a low mean?</strong> </p>

<hr>

<p>In response to @Jarod42's comment who said</p>

<blockquote>
  <p>Wonder if using fewer bits don't break equiprobability...</p>
</blockquote>

<p>I don't think it would break equiprobability. In a vague attempt to test it, I consider the same question with a simple bernoulli distribution. I am sampling true with a probability <code>1/2^4</code> and sampling false with a probability <code>1 - 1/2^4</code>. The function <code>drawWithoutWastingBits</code> stops as soon as it sees a true in the cache and the function <code>drawWastingBits</code> consumes 4 bits whatever these bits are.</p>

<pre><code>#include &lt;iostream&gt;
#include &lt;vector&gt;
#include &lt;string&gt;
#include &lt;algorithm&gt;
#include &lt;random&gt;

bool drawWithoutWastingBits(std::vector&lt;bool&gt;&amp; cache, size_t&amp; cache_index)
{
    /* 
        Get a true with probability 1/2^4 (=1/16=0.0625) and a false otherwise
    */

    size_t nbTrues = 0;
    while (cache[cache_index])
    {
        ++nbTrues;
        ++cache_index;
        if (nbTrues == 4)
        {
            return true;
        }
    }
    ++cache_index;
    return false;
}


bool drawWastingBits(std::vector&lt;bool&gt;&amp; cache, size_t&amp; cache_index)
{
    /* 
        Get a true with probability 1/2^4 (=1/16=0.0625) and a false otherwise
    */

    bool isAnyTrue = false;
    for (size_t i = 0 ; i &lt; 4; ++i)
    {
        if (cache[cache_index])
        {
            isAnyTrue = true;
        }
        ++cache_index;
    }
    return !isAnyTrue;
}

int main()
{
    /*
        Just cache a lot of bits in advance in `cache`. The same sequence of bits will be used by both function.
        I am just caching way enough bits to make sure they don't run out of bits below
        I made sure to have the same number of zeros and ones so that any deviation is caused by the methodology and not by the RNG
    */

    // Produce cache
    std::vector&lt;bool&gt; cache;
    size_t nbBitsToCache = 1e7;
    cache.reserve(nbBitsToCache);
    for (size_t i = 0 ; i &lt; nbBitsToCache/2 ; ++i)
    {
        cache.push_back(false);
        cache.push_back(true);
    }
    // Shuffle cache
    {
        std::mt19937 mt(std::random_device{}());
        std::shuffle(cache.begin(), cache.end(), mt);
    }


    // Draw without wasting bits
    {
        size_t nbDraws = 1e6;
        size_t cache_index = 0;
        std::pair&lt;size_t, size_t&gt; outcomes = {0,0};
        for (size_t r = 0 ; r &lt; nbDraws ; ++r)
        {
            drawWithoutWastingBits(cache, cache_index) ? ++outcomes.first : ++outcomes.second;
            assert(cache_index &lt;= cache.size());
        }   

        assert(outcomes.first + outcomes.second == nbDraws);
        std::cout &lt;&lt; "Draw Without Wasting Bits: prob true = " &lt;&lt; (double)outcomes.first / nbDraws &lt;&lt; "\n";
    }


    // Draw wasting bits
    {
        size_t nbDraws = 1e6;
        size_t cache_index = 0;
        std::pair&lt;size_t, size_t&gt; outcomes = {0,0};
        for (size_t r = 0 ; r &lt; nbDraws ; ++r)
        {
            drawWastingBits(cache, cache_index) ? ++outcomes.first : ++outcomes.second;
            assert(cache_index &lt;= cache.size());
        }   

        assert(outcomes.first + outcomes.second == nbDraws);
        std::cout &lt;&lt; "Draw Wit Wasting Bits: prob true = " &lt;&lt; (double)outcomes.first / nbDraws &lt;&lt; "\n";
    }
}
</code></pre>

<p>Possible output</p>

<pre><code>Draw Without Wasting Bits: prob true = 0.062832
Draw Wit Wasting Bits: prob true = 0.062363
</code></pre>