src/DIGRand/DIGPoisson.cpp

   1 /*
   2  GENERATE POISSON RANDOM NUMBER IF alambda IS LESS THAN xover
   3  OR ELSE GENERATE GAUSSIAN RANDOM NUMBER
   4 */
   5
   6 #include <cstdio>
   7 #include <cmath>
   8
   9 //#include "DIGRand.h"
  10 #include <DIGRand/DIGRand.h>
  11 //#include "DIGGauss.h"
  12 #include <DIGRand/DIGGauss.h>
  13 //#include "DIGPoisson.h"
  14 #include <DIGRand/DIGPoisson.h>
  15
  16
  17 int Poisson(double pALambda)
  18 {
  19   static const double xover = 36.0;
  20
  21   int n;
  22
  23   double a, e, y;
  24
  25
  26   if(pALambda >= xover) {
  27     int g = (int) (Gauss(pALambda, (double) sqrt(pALambda)) + 0.5);
  28     n = (g > 0) ? g: 0;
  29   }
  30   else {
  31     y = Rand();
  32     n = 0;
  33     e = (double) exp(-pALambda);
  34     a = y;
  35
  36     while(a >= e) {
  37       n++;
  38       y = Rand();
  39       a *= y;
  40     }
  41   }
  42   return n;
  43 }