diff options
author | Nicholas <nbnoll@eml.cc> | 2021-11-12 09:22:01 -0800 |
---|---|---|
committer | Nicholas <nbnoll@eml.cc> | 2021-11-12 09:22:01 -0800 |
commit | ce05175372a9ddca1a225db0765ace1127a39293 (patch) | |
tree | 5988b4d4f6b402e4953945886fc90aae11203df6 /src/base/rng/poisson.c | |
parent | b375f3cdedb5b0e08745d100b40e38d2f8396a58 (diff) |
chore: simplified organizational structurelaptop
Diffstat (limited to 'src/base/rng/poisson.c')
-rw-r--r-- | src/base/rng/poisson.c | 126 |
1 files changed, 126 insertions, 0 deletions
diff --git a/src/base/rng/poisson.c b/src/base/rng/poisson.c new file mode 100644 index 0000000..3ec15c9 --- /dev/null +++ b/src/base/rng/poisson.c @@ -0,0 +1,126 @@ +#include "internal.h" + +/* + * Ahrens, J. H., & Dieter, U. (1982). + * Computer Generation of Poisson Deviates from Modified Normal Distributions. + */ +static double factorial[10] = {1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880.}; +static double coeffs[9] = { + -.500000000, +.333333333, -.249999856, + +.200011780, -.166684875, +.142187833, + -.124196313, +.125005956, -.114265030, +}; + +static inline +double +log1pmx(double x, double off) +{ + int i; + double r, t; + + if(-0.25 < x && x < 0.25) { + r = 0; + t = 1; + for(i=0;i<arrlen(coeffs);i++) { + r += coeffs[i]*t; + t *= x; + } + + return x*x*r; + } + return log(1+x) - off; +} + +static inline +double +procf(double mu, double s, int64 K, double *px, double *py, double *fx, double *fy) +{ + double d, V, X; + double w, b1, b2, c1, c2, c3, c0, c; + + w = 0.3989422804014327/s; + b1 = 0.041666666666666664/mu; + b2 = 0.3*b1*b1; + c3 = 0.14285714285714285*b1*b2; + c2 = b2 - 15.*c3; + c1 = b1 - 6.*b2 + 45.*c3; + c0 = 1 - b1 + 3.*b2 - 15.*c3; + c = .1069/mu; + + if(K < 10) { + *px = -mu; + *py = pow(mu,K) / factorial[K]; + }else{ + d = 0.08333333333333333/K; + d = d - 4.8*d*d*d; + V = (mu - K) / K; + + *px = K*log1pmx(V,mu-K) - d; + *py = 0.3989422804014327/sqrt(K); + } + + X = (K - mu + 0.5)/s; + *fx = -0.5*X*X; + *fy = w*(((c3*X*X + c2)*X*X + c1)*X*X + c0); + + return c; +} + +static inline +uint64 +bigpoisson(double mu) +{ + int64 L,K; + double G,s,d,U,E,T; + double px,py,fx,fy,c; + + s = sqrt(mu); + d = 6*mu*mu; + L = floor(mu - 1.1484); + +stepN: + G = mu + s*rng·normal(); + K = floor(G); + if(K<0) + goto stepP; +stepI: + if(K>=L) + return K; +stepS: + U = rng·random(); + if(d*U >= (mu-K)*(mu-K)*(mu-K)) + return K; +stepP: + if(G < 0) + goto stepE; +stepQ: + c = procf(mu, s, K, &px, &py, &fx, &fy); +stepE: + E = rng·exponential(1.0); + U = rng·random(); + U = U + U - 1; + T = 1.8 + copysign(E,U); + if(T < 0.6744) + goto stepE; + K = floor(mu + s*T); + c = procf(mu, s, K, &px, &py, &fx, &fy); +stepH: + if(c*fabs(U) > (py*exp(px + E) - fy*exp(fx + E))) + goto stepE; + return K; +} + +uint64 +rng·poisson(double mean) +{ + int64 n; + double z; + + if(mean<10.0) { + for(n=0, z=rng·exponential(1.0); z<mean; ++n, z+=rng·exponential(1.0)) + ; + return n; + } + + return bigpoisson(mean); +} |