1 /*
2 * Copyright (C) 1994. James Darrell McCauley. (darrell@mccauley-usa.com)
3 * http://mccauley-usa.com/
4 *
5 * This program is free software under the GPL (>=v2)
6 * Read the file GPL.TXT coming with GRASS for details.
7 */
8
9 #include <stdio.h>
10 #include <math.h>
11 #include "zufall.h"
12
13 /*-Poisson generator for distribution function of p's:
14 * q(mu,p) = exp(-mu) mu**p/p!
15 */
16
fische(int n,double * mu,int * p)17 int fische(int n, double *mu, int *p)
18 {
19 int left, indx[1024], i, k, nsegs, p0, ii, jj, nl0;
20 double q[1024], u[1024], q0, pmu;
21
22 --p;
23
24 if (n <= 0)
25 return 0;
26
27 pmu = exp(-(*mu));
28 p0 = 0;
29
30 nsegs = (n - 1) / 1024;
31 left = n - (nsegs << 10);
32 ++nsegs;
33 nl0 = left;
34
35 for (k = 1; k <= nsegs; ++k) {
36 for (i = 1; i <= left; ++i) {
37 indx[i - 1] = i;
38 p[p0 + i] = 0;
39 q[i - 1] = 1.0;
40 }
41
42 do { /* Begin iterative loop on segment of p's */
43 zufall(left, u); /* Get the needed uniforms */
44 jj = 0;
45
46 for (i = 1; i <= left; ++i) {
47 ii = indx[i - 1];
48 q0 = q[ii - 1] * u[i - 1];
49 q[ii - 1] = q0;
50 if (q0 > pmu) {
51 indx[jj++] = ii;
52 ++p[p0 + ii];
53 }
54 }
55 left = jj; /* any left in this segment? */
56 }
57 while (left > 0);
58
59 p0 += nl0;
60 nl0 = left = 1024;
61 }
62 return 0;
63 } /* fische */
64