/********************************************************************/ /* rand.c */ /* */ /* Routine to return a random number from a Poisson distribution */ /* Includes the following routines from Numerical Recipes in C, */ /* p 222-223, by Press WH, et al. Cambridge University Press (C) */ /* 1988: */ /* */ /* expdev(), poidev() and gammln() */ /* */ /* Please ensure you have the necessary licenses to use these. */ /* */ /********************************************************************/ #include #include #include #include #define PI 3.141592654 double gammln(double x); double poidev(double xm, long *idum); double expdev (double xm, long *idum); double expdev(double xm, long *idum) { double dum; do dum = ran2(idum); while (dum == 0.0); return (-xm*log(dum)); } double poidev(double xm, long *idum) { static double sq, alxm, g, oldm=(-1.0); double em, t, y; if (xm < 12.0) { if (xm != oldm) { oldm = xm; g = exp(-xm); } em = -1; t = 1.0; do { em += 1.0; t *= ran2(idum); } while (t > g); } else { if (xm != oldm) { oldm = xm; sq = sqrt(2.0*xm); alxm = log(xm); g = xm * alxm - gammln(xm+1.0); } do { do { y = tan(PI*ran2(idum)); em = sq * y + xm; } while (em < 0.0); em = floor(em); t = 0.9 * (1.0 + y*y) * exp(em * alxm - gammln(em+1.0) - g); } while (ran2(idum) > t); } return em; } double gammln(double xx) { double x, tmp, ser; static double cof[6]={76.18009173,-86.50532033,24.01409822, -1.231739516,0.120858003e-2,-0.536382e-5}; int j; x = xx - 1.0; tmp = x + 5.5; tmp -= (x + 0.5) * log(tmp); ser = 1.0; for (j=0; j<=5; j++) { x += 1.0; ser += cof[j]/x; } return -tmp + log(2.50662827465 * ser); }