1.1 --- a/lemon/random.h	Mon Mar 17 18:31:52 2008 +0000
     1.2 +++ b/lemon/random.h	Mon Mar 17 19:21:27 2008 +0000
     1.3 @@ -803,6 +803,29 @@
     1.4        return exponential(gamma(k,1.0/x_min));
     1.5      }  
     1.6        
     1.7 +    /// Poisson distribution
     1.8 +
     1.9 +    /// This function generates a Poisson distribution random number with
    1.10 +    /// parameter \c lambda.
    1.11 +    /// 
    1.12 +    /// The probability mass function of this distribusion is
    1.13 +    /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
    1.14 +    /// \note The algorithm is taken from the book of Donald E. Knuth titled
    1.15 +    /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
    1.16 +    /// return value.
    1.17 +    
    1.18 +    int poisson(double lambda)
    1.19 +    {
    1.20 +      const double l = std::exp(-lambda);
    1.21 +      int k=0;
    1.22 +      double p = 1.0;
    1.23 +      do {
    1.24 +	k++;
    1.25 +	p*=real<double>();
    1.26 +      } while (p>=l);
    1.27 +      return k-1;
    1.28 +    }  
    1.29 +      
    1.30      ///@}
    1.31      
    1.32      ///\name Two dimensional distributions

     2.1 --- a/test/random_test.cc	Mon Mar 17 18:31:52 2008 +0000
     2.2 +++ b/test/random_test.cc	Mon Mar 17 19:21:27 2008 +0000
     2.3 @@ -33,4 +33,5 @@
     2.4    a=lemon::rnd.gamma(4);
     2.5    //Does gamma work with integer k?
     2.6    a=lemon::rnd.gamma(4.0,0);
     2.7 +  a=lemon::rnd.poisson(.5);
     2.8  }