/* ranpoi.c: random # ~ Poisson(lambda) */
/***********************************************************/
#include <math.h>
#include "ran.h"

int  RanPoi(double lambda)
/* returns Poisson(lambda) */
{
  int y = 0;
  double s = 0;

  if (lambda < 0) return -1;

  else if (lambda < 16) {       /* small lambda */
    while ((s += RanExp()) < lambda) y++;
    return y;
  }
  else {                        /* large lambda */
    y = 0.5 + lambda + sqrt(lambda)*RanNor(); /* rounded */
    return (y<0 ? 0: y);
  }
}