; Code developed by Melanie Mitchell and Stephanie Forrest for their Royal Road
; GA work.

;------------------------------------------------------------------------

(defun poisson (mean &optional (compute-lambda-term? t))
; Returns a random number (an integer) from a poisson distribution.
; The mean of the distribution is passed into the function as well
; as an optional flag (compute-lambda-term?).  If the flag is t,
; then *lambda-term* is set to expt(- mean).   Otherwise,
; *lambda-term* is assumed to have the value exp (- mean).  Use of
; this term avoids unnecessary recalculation of the exponent.  Note,
; *lambda-term* is not identical with lambda of the poisson
; distribution, since lambda = mean.

    (let ((uniform (new-random-float))
          term
          (k-value 0)
          cum)
      (declare (special *lambda-term*))
        (if (< mean 0) (error "Error in Poisson!!: bad mean"))

        (if compute-lambda-term?
          (setf *lambda-term* (exp (- mean))))

        ;; Compute the cumulative density function (c.d.f.) for successive
        ;; values of k-value.  When the c.d.f. becomes greater than the
        ;; random variable chosen from the Uniform distribution (uniforms),
        ;; we quit.  Since k must be non-negative, and since
        ;; prob(k=0) = lambda_term, k will start out at 0 with the other values
        ;; initialized accordingly.

        (setf term *lambda-term*)
        (setf cum *lambda-term*)
        (loop
            (if (or (> cum uniform) (<= term 0.0000005)) ; Cutoff for high end
                                                         ; of distribution.
                (return k-value))

           (incf k-value)
           (setf term (/ (* term  mean) k-value))
           (setf cum (+ cum term)))
    k-value))

;------------------------------------------------------------------------

This is written in common lisp.  Steele's book called
Common Lisp is what I had been using for reference.
new-random-float is a function that generates a random float.