covid-sim
Loading...
Searching...
No Matches
Rand.cpp
1#include <stdio.h>
2#include <cmath>
3#include <algorithm>
4#include "Rand.h"
5#include "MachineDefines.h"
6#include "Constants.h"
7#include "Error.h"
8#ifdef _OPENMP
9#include <omp.h>
10#endif
11
12/* RANDLIB static variables */
13int32_t* Xcg1, *Xcg2;
14int **SamplingQueue = nullptr;
15
19
20double ranf(void)
21{
22 int32_t k, s1, s2, z;
23 unsigned int curntg;
24
25#ifdef _OPENMP
26 curntg = CACHE_LINE_SIZE * omp_get_thread_num();
27#else
28 curntg = 0;
29#endif
30 s1 = Xcg1[curntg];
31 s2 = Xcg2[curntg];
32 k = s1 / 53668;
33 s1 = Xa1 * (s1 - k * 53668) - k * 12211;
34 if (s1 < 0) s1 += Xm1;
35 k = s2 / 52774;
36 s2 = Xa2 * (s2 - k * 52774) - k * 3791;
37 if (s2 < 0) s2 += Xm2;
38 Xcg1[curntg] = s1;
39 Xcg2[curntg] = s2;
40 z = s1 - s2;
41 if (z < 1) z += (Xm1 - 1);
42 return ((double)z) / Xm1;
43}
44
45double ranf_mt(int tn)
46{
47 int32_t k, s1, s2, z;
48 int curntg;
49
50 curntg = CACHE_LINE_SIZE * tn;
51 s1 = Xcg1[curntg];
52 s2 = Xcg2[curntg];
53 k = s1 / 53668;
54 s1 = Xa1 * (s1 - k * 53668) - k * 12211;
55 if (s1 < 0) s1 += Xm1;
56 k = s2 / 52774;
57 s2 = Xa2 * (s2 - k * 52774) - k * 3791;
58 if (s2 < 0) s2 += Xm2;
59 Xcg1[curntg] = s1;
60 Xcg2[curntg] = s2;
61 z = s1 - s2;
62 if (z < 1) z += (Xm1 - 1);
63 return ((double)z) / Xm1;
64}
65
66void setall(int32_t *pseed1, int32_t *pseed2)
67/*
68**********************************************************************
69 void setall(int32_t iseed1,int32_t iseed2)
70 SET ALL random number generators
71 Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
72 initial seeds of the other generators are set accordingly, and
73 all generators states are set to these seeds.
74 This is a transcription from Pascal to C++ of routine
75 Set_Initial_Seed from the paper
76 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
77 with Splitting Facilities." ACM Transactions on Mathematical
78 Software, 17:98-111 (1991)
79 https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.149.9439
80 Arguments
81 iseed1 -> First of two integer seeds
82 iseed2 -> Second of two integer seeds
83**********************************************************************
84*/
85{
86 int g;
87
88 int32_t iseed1 = *pseed1;
89 int32_t iseed2 = *pseed2;
90
91 for (g = 0; g < MAX_NUM_THREADS; g++) {
92 *(Xcg1 + g * CACHE_LINE_SIZE) = iseed1 = mltmod(Xa1vw, iseed1, Xm1);
93 *(Xcg2 + g * CACHE_LINE_SIZE) = iseed2 = mltmod(Xa2vw, iseed2, Xm2);
94 }
95
96 *pseed1 = iseed1;
97 *pseed2 = iseed2;
98}
99
100int32_t mltmod(int32_t a, int32_t s, int32_t m)
101/*
102**********************************************************************
103 int32_t mltmod(int32_t a, int32_t s, int32_t m)
104 Returns (a * s) MOD m
105 This is a transcription from Pascal to C++ of routine
106 MULtMod_Decompos from the paper
107 L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
108 with Splitting Facilities." ACM Transactions on Mathematical
109 Software, 17:98-111 (1991)
110 https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.149.9439
111**********************************************************************
112*/
113{
114 const int32_t h = 32768;
115 int32_t a0, a1, k, p, q, qh, rh;
116 /*
117 H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
118 machine. On a different machine recompute H
119 */
120 if (a <= 0 || a >= m || s <= 0 || s >= m) {
121 fputs(" a, m, s out of order in mltmod - ABORT!\n", stderr);
122 fprintf(stderr, " a = %12d s = %12d m = %12d\n", a, s, m);
123 fputs(" mltmod requires: 0 < a < m; 0 < s < m\n", stderr);
124 exit(1);
125 }
126
127 if (a < h) {
128 a0 = a;
129 p = 0;
130 } else {
131 a1 = a / h;
132 a0 = a - h * a1;
133 qh = m / h;
134 rh = m - h * qh;
135 if (a1 >= h) { // a2 == 1
136 a1 -= h;
137 k = s / qh;
138 p = h * (s - k * qh) - k * rh;
139 while (p < 0) {
140 p += m;
141 }
142 } else {
143 p = 0;
144 }
145 // p == (a2 * s * h) MOD m
146 if (a1 != 0) {
147 q = m / a1;
148 k = s / q;
149 p -= (k * (m - a1 * q));
150 if (p > 0) p -= m;
151 p += (a1 * (s - k * q));
152 while (p < 0) {
153 p += m;
154 }
155 }
156 // p == ((a2 * h + a1) * s) MOD m
157 k = p / qh;
158 p = h * (p - k * qh) - k * rh;
159 while (p < 0) {
160 p += m;
161 }
162 }
163 // p == ((a2 * h + a1) * h * s) MOD m
164 if (a0 != 0) {
165 q = m / a0;
166 k = s / q;
167 p -= (k * (m - a0 * q));
168 if (p > 0) p -= m;
169 p += (a0 * (s - k * q));
170 while (p < 0) {
171 p += m;
172 }
173 }
174 return p;
175}
176
177int32_t ignbin(int32_t n, double pp)
178{
179 /*
180**********************************************************************
181 int32_t ignbin(int32_t n,double pp)
182 GENerate BINomial random deviate
183 Function
184 Generates a single random deviate from a binomial
185 distribution whose number of trials is N and whose
186 probability of an event in each trial is P.
187 Arguments
188 n --> The number of trials in the binomial distribution
189 from which a random deviate is to be generated.
190 JJV (N >= 0)
191 pp --> The probability of an event in each trial of the
192 binomial distribution from which a random deviate
193 is to be generated.
194 JJV (0.0 <= PP <= 1.0)
195 ignbin <-- A random deviate yielding the number of events
196 from N independent trials, each of which has
197 a probability of event P.
198 Method
199 This is algorithm BTPE from:
200 Kachitvichyanukul, V. and Schmeiser, B. W.
201 Binomial Random Variate Generation.
202 Communications of the ACM, 31, 2
203 (February, 1988) 216.
204**********************************************************************
205 SUBROUTINE BTPEC(N,PP,ISEED,JX)
206 BINOMIAL RANDOM VARIATE GENERATOR
207 MEAN .LT. 30 -- INVERSE CDF
208 MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA
209 FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE
210 (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
211 THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
212 BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
213 BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE
214 RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
215 USABLE ALGORITHM.
216 REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
217 "BINOMIAL RANDOM VARIATE GENERATION,"
218 COMMUNICATIONS OF THE ACM, FORTHCOMING
219 WRITTEN: SEPTEMBER 1980.
220 LAST REVISED: MAY 1985, JULY 1987
221 REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER
222 GENERATOR
223 ARGUMENTS
224 N : NUMBER OF BERNOULLI TRIALS (INPUT)
225 PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
226 ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
227 JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)
228 VARIABLES
229 PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
230 NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
231 XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
232 P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
233 FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
234 M: INTEGER VALUE OF THE CURRENT MODE
235 FM: FLOATING POINT VALUE OF THE CURRENT MODE
236 XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
237 P1: AREA OF THE TRIANGLE
238 C: HEIGHT OF THE PARALLELOGRAMS
239 XM: CENTER OF THE TRIANGLE
240 XL: LEFT END OF THE TRIANGLE
241 XR: RIGHT END OF THE TRIANGLE
242 AL: TEMPORARY VARIABLE
243 XLL: RATE FOR THE LEFT EXPONENTIAL TAIL
244 XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL
245 P2: AREA OF THE PARALLELOGRAMS
246 P3: AREA OF THE LEFT EXPONENTIAL TAIL
247 P4: AREA OF THE RIGHT EXPONENTIAL TAIL
248 U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
249 FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
250 FROM THE REGION
251 V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
252 (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
253 REJECT THE CANDIDATE VALUE
254 IX: INTEGER CANDIDATE VALUE
255 X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
256 AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
257 K: ABSOLUTE VALUE OF (IX-M)
258 F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
259 ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
260 ALSO USED IN THE INVERSE TRANSFORMATION
261 R: THE RATIO P/Q
262 G: CONSTANT USED IN CALCULATION OF PROBABILITY
263 MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
264 OF F WHEN IX IS GREATER THAN M
265 IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
266 CALCULATION OF F WHEN IX IS LESS THAN M
267 I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
268 AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
269 YNORM: LOGARITHM OF NORMAL BOUND
270 ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
271 X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
272 USED IN THE FINAL ACCEPT/REJECT TEST
273 QN: PROBABILITY OF NO SUCCESS IN N TRIALS
274 REMARK
275 IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
276 SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME
277 COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
278 ARE NOT INVOLVED.
279 ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
280 GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
281 TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
282**********************************************************************
283*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
284*/
285/* JJV changed initial values to ridiculous values */
286 double psave = -1.0E37;
287 int32_t nsave = -214748365;
288 int32_t ignbin, i, ix, ix1, k, m, mp, T1;
289 double al, alv, amaxp, c, f, f1, f2, ffm, fm, g, p, p1, p2, p3, p4, q, qn, r, u, v, w, w2, x, x1,
290 x2, xl, xll, xlr, xm, xnp, xnpq, xr, ynorm, z, z2;
291
292 /*
293 *****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
294 JJV added checks to ensure 0.0 <= PP <= 1.0
295 */
296 if (pp < 0.0F) ERR_CRITICAL("PP < 0.0 in IGNBIN");
297 if (pp > 1.0F) ERR_CRITICAL("PP > 1.0 in IGNBIN");
298 psave = pp;
299 p = std::min(psave, 1.0 - psave);
300 q = 1.0 - p;
301
302 /*
303 JJV added check to ensure N >= 0
304 */
305 if (n < 0L) ERR_CRITICAL("N < 0 in IGNBIN");
306 xnp = n * p;
307 nsave = n;
308 if (xnp < 30.0) goto S140;
309 ffm = xnp + p;
310 m = (int32_t)ffm;
311 fm = m;
312 xnpq = xnp * q;
313 p1 = (int32_t)(2.195 * sqrt(xnpq) - 4.6 * q) + 0.5;
314 xm = fm + 0.5;
315 xl = xm - p1;
316 xr = xm + p1;
317 c = 0.134 + 20.5 / (15.3 + fm);
318 al = (ffm - xl) / (ffm - xl * p);
319 xll = al * (1.0 + 0.5 * al);
320 al = (xr - ffm) / (xr * q);
321 xlr = al * (1.0 + 0.5 * al);
322 p2 = p1 * (1.0 + c + c);
323 p3 = p2 + c / xll;
324 p4 = p3 + c / xlr;
325S30:
326 /*
327 *****GENERATE VARIATE
328 */
329 u = ranf() * p4;
330 v = ranf();
331 /*
332 TRIANGULAR REGION
333 */
334 if (u > p1) goto S40;
335 ix = (int32_t)(xm - p1 * v + u);
336 goto S170;
337S40:
338 /*
339 PARALLELOGRAM REGION
340 */
341 if (u > p2) goto S50;
342 x = xl + (u - p1) / c;
343 v = v * c + 1.0 - std::abs(xm - x) / p1;
344 if (v > 1.0 || v <= 0.0) goto S30;
345 ix = (int32_t)x;
346 goto S70;
347S50:
348 /*
349 LEFT TAIL
350 */
351 if (u > p3) goto S60;
352 ix = (int32_t)(xl + log(v) / xll);
353 if (ix < 0) goto S30;
354 v *= ((u - p2) * xll);
355 goto S70;
356S60:
357 /*
358 RIGHT TAIL
359 */
360 ix = (int32_t)(xr - log(v) / xlr);
361 if (ix > n) goto S30;
362 v *= ((u - p3) * xlr);
363S70:
364 /*
365 *****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
366 */
367 k = std::abs(ix - m);
368 if (k > 20 && k < xnpq / 2 - 1) goto S130;
369 /*
370 EXPLICIT EVALUATION
371 */
372 f = 1.0;
373 r = p / q;
374 g = (n + 1) * r;
375 T1 = m - ix;
376 if (T1 < 0) goto S80;
377 else if (T1 == 0) goto S120;
378 else goto S100;
379S80:
380 mp = m + 1;
381 for (i = mp; i <= ix; i++) f *= (g / i - r);
382 goto S120;
383S100:
384 ix1 = ix + 1;
385 for (i = ix1; i <= m; i++) f /= (g / i - r);
386S120:
387 if (v <= f) goto S170;
388 goto S30;
389S130:
390 /*
391 SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
392 */
393 amaxp = k / xnpq * ((k * (k / 3.0 + 0.625) + 0.1666666666666) / xnpq + 0.5);
394 ynorm = -(k * k / (2.0 * xnpq));
395 alv = log(v);
396 if (alv < ynorm - amaxp) goto S170;
397 if (alv > ynorm + amaxp) goto S30;
398 /*
399 STIRLING'S FORMULA TO MACHINE ACCURACY FOR
400 THE FINAL ACCEPTANCE/REJECTION TEST
401 */
402 x1 = ix + 1.0;
403 f1 = fm + 1.0;
404 z = n + 1.0 - fm;
405 w = n - ix + 1.0;
406 z2 = z * z;
407 x2 = x1 * x1;
408 f2 = f1 * f1;
409 w2 = w * w;
410 if (alv <= xm * log(f1 / x1) + (n - m + 0.5) * log(z / w) + (ix - m) * log(w * p / (x1 * q)) + (13860.0 -
411 (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 - (462.0 -
412 (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 - (462.0 - (132.0 -
413 (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0
414 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.0) goto S170;
415 goto S30;
416S140:
417 /*
418 INVERSE CDF LOGIC FOR MEAN LESS THAN 30
419 */
420 qn = pow(q, (double)n);
421 r = p / q;
422 g = r * (n + 1);
423S150:
424 ix = 0;
425 f = qn;
426 u = ranf();
427S160:
428 if (u < f) goto S170;
429 if (ix > 110) goto S150;
430 u -= f;
431 ix += 1;
432 f *= (g / ix - r);
433 goto S160;
434S170:
435 if (psave > 0.5) ix = n - ix;
436 ignbin = ix;
437 return ignbin;
438}
439
440int32_t ignpoi(double mu)
441/*
442**********************************************************************
443 int32_t ignpoi(double mu)
444 GENerate POIsson random deviate
445 Function
446 Generates a single random deviate from a Poisson
447 distribution with mean MU.
448 Arguments
449 mu --> The mean of the Poisson distribution from which
450 a random deviate is to be generated.
451 (mu >= 0.0)
452 ignpoi <-- The random deviate.
453 Method
454 Renames KPOIS from TOMS as slightly modified by BWB to use RANF
455 instead of SUNIF.
456 For details see:
457 Ahrens, J.H. and Dieter, U.
458 Computer Generation of Poisson Deviates
459 From Modified Normal Distributions.
460 ACM Trans. Math. Software, 8, 2
461 (June 1982),163-179
462**********************************************************************
463**********************************************************************
464
465
466 P O I S S O N DISTRIBUTION
467
468
469**********************************************************************
470**********************************************************************
471
472 FOR DETAILS SEE:
473
474 AHRENS, J.H. AND DIETER, U.
475 COMPUTER GENERATION OF POISSON DEVIATES
476 FROM MODIFIED NORMAL DISTRIBUTIONS.
477 ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
478
479 (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
480
481**********************************************************************
482 INTEGER FUNCTION IGNPOI(IR,MU)
483 INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
484 MU=MEAN MU OF THE POISSON DISTRIBUTION
485 OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
486 MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
487 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
488 COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
489 SEPARATION OF CASES A AND B
490*/
491{
492 extern double fsign(double num, double sign);
493 static double a0 = -0.5;
494 static double a1 = 0.3333333;
495 static double a2 = -0.2500068;
496 static double a3 = 0.2000118;
497 static double a4 = -0.1661269;
498 static double a5 = 0.1421878;
499 static double a6 = -0.1384794;
500 static double a7 = 0.125006;
501 /* JJV changed the initial values of MUPREV and MUOLD */
502 static double fact[10] = {
503 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0
504 };
505 /* JJV added ll to the list, for Case A */
506 int32_t ignpoi, j, k, kflag, l, ll, m;
507 double b1, b2, c, c0, c1, c2, c3, d, del, difmuk, e, fk, fx, fy, g, omega, p, p0, px, py, q, s, t, u, v, x, xx, pp[35];
508
509 if (mu < 10.0) goto S120;
510 /*
511 C A S E A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED)
512 JJV changed l in Case A to ll
513 */
514 s = sqrt(mu);
515 d = 6.0 * mu * mu;
516 /*
517 THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
518 PROBABILITIES FK WHENEVER K >= M(MU). LL=IFIX(MU-1.1484)
519 IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
520 */
521 ll = (int32_t)(mu - 1.1484);
522
523 /*
524 STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
525 */
526 g = mu + s * snorm();
527 if (g < 0.0) goto S20;
528 ignpoi = (int32_t)(g);
529 /*
530 STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
531 */
532 if (ignpoi >= ll) return ignpoi;
533 /*
534 STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
535 */
536 fk = (double)ignpoi;
537 difmuk = mu - fk;
538 u = ranf();
539 if (d * u >= difmuk * difmuk * difmuk) return ignpoi;
540S20:
541 /*
542 STEP P. PREPARATIONS FOR STEPS Q AND H.
543 (RECALCULATIONS OF PARAMETERS IF NECESSARY)
544 .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
545 THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
546 APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
547 C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
548 */
549 omega = 0.3989423 / s;
550 b1 = 4.166667E-2 / mu;
551 b2 = 0.3 * b1 * b1;
552 c3 = 0.1428571 * b1 * b2;
553 c2 = b2 - 15.0 * c3;
554 c1 = b1 - 6.0 * b2 + 45.0 * c3;
555 c0 = 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
556 c = 0.1069 / mu;
557
558 if (g < 0.0) goto S50;
559 /*
560 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
561 */
562 kflag = 0;
563 goto S70;
564S40:
565 /*
566 STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
567 */
568 if (fy - u * fy <= py * exp(px - fx)) return ignpoi;
569S50:
570 /*
571 STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
572 DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
573 (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
574 */
575 e = sexpo();
576 u = ranf();
577 u += (u - 1.0);
578 t = 1.8 + fsign(e, u);
579 if (t <= -0.6744) goto S50;
580 ignpoi = (int32_t)(mu + s * t);
581 fk = (double)ignpoi;
582 difmuk = mu - fk;
583 /*
584 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
585 */
586 kflag = 1;
587 goto S70;
588S60:
589 /*
590 STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
591 */
592 if (c * fabs(u) > py * exp(px + e) - fy * exp(fx + e)) goto S50;
593 return ignpoi;
594S70:
595 /*
596 STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
597 CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
598 */
599 if (ignpoi >= 10) goto S80;
600 px = -mu;
601 py = pow(mu, (double)ignpoi) / *(fact + ignpoi);
602 goto S110;
603S80:
604 /*
605 CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
606 A0-A7 FOR ACCURACY WHEN ADVISABLE
607 .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
608 */
609 del = 8.333333E-2 / fk;
610 del -= (4.8 * del * del * del);
611 v = difmuk / fk;
612 if (fabs(v) <= 0.25) goto S90;
613 px = fk * log(1.0 + v) - difmuk - del;
614 goto S100;
615S90:
616 px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v + a0) - del;
617S100:
618 py = 0.3989423 / sqrt(fk);
619S110:
620 x = (0.5 - difmuk) / s;
621 xx = x * x;
622 fx = -0.5 * xx;
623 fy = omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
624 if (kflag <= 0) goto S40;
625 goto S60;
626S120:
627 /*
628 C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
629 JJV changed MUPREV assignment to initial value
630 */
631 m = std::max(INT32_C(1), (int32_t)(mu));
632 l = 0;
633 p = exp(-mu);
634 q = p0 = p;
635S130:
636 /*
637 STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
638 */
639 u = ranf();
640 ignpoi = 0;
641 if (u <= p0) return ignpoi;
642 /*
643 STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
644 PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
645 (0.458=PP(9) FOR MU=10)
646 */
647 if (l == 0) goto S150;
648 j = 1;
649 if (u > 0.458) j = std::min(l, m);
650 for (k = j; k <= l; k++) {
651 if (u <= *(pp + k - 1)) goto S180;
652 }
653 if (l == 35) goto S130;
654S150:
655 /*
656 STEP C. CREATION OF NEW POISSON PROBABILITIES P
657 AND THEIR CUMULATIVES Q=PP(K)
658 */
659 l += 1;
660 for (k = l; k <= 35; k++) {
661 p = p * mu / (double)k;
662 q += p;
663 *(pp + k - 1) = q;
664 if (u <= q) goto S170;
665 }
666 l = 35;
667 goto S130;
668S170:
669 l = k;
670S180:
671 ignpoi = k;
672 return ignpoi;
673}
674
675int32_t ignpoi_mt(double mu, int tn)
676/*
677**********************************************************************
678int32_t ignpoi_mt(double mu)
679GENerate POIsson random deviate
680Function
681Generates a single random deviate from a Poisson
682distribution with mean MU.
683Arguments
684mu --> The mean of the Poisson distribution from which
685a random deviate is to be generated.
686(mu >= 0.0)
687ignpoi_mt <-- The random deviate.
688Method
689Renames KPOIS from TOMS as slightly modified by BWB to use RANF
690instead of SUNIF.
691For details see:
692Ahrens, J.H. and Dieter, U.
693Computer Generation of Poisson Deviates
694From Modified Normal Distributions.
695ACM Trans. Math. Software, 8, 2
696(June 1982),163-179
697**********************************************************************
698**********************************************************************
699
700
701P O I S S O N DISTRIBUTION
702
703
704**********************************************************************
705**********************************************************************
706
707FOR DETAILS SEE:
708
709AHRENS, J.H. AND DIETER, U.
710COMPUTER GENERATION OF POISSON DEVIATES
711FROM MODIFIED NORMAL DISTRIBUTIONS.
712ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
713
714(SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
715
716**********************************************************************
717INTEGER FUNCTION IGNPOI(IR,MU)
718INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
719MU=MEAN MU OF THE POISSON DISTRIBUTION
720OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
721MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
722TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
723COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
724SEPARATION OF CASES A AND B
725*/
726{
727 extern double fsign(double num, double sign);
728 double a0 = -0.5;
729 double a1 = 0.3333333;
730 double a2 = -0.2500068;
731 double a3 = 0.2000118;
732 double a4 = -0.1661269;
733 double a5 = 0.1421878;
734 double a6 = -0.1384794;
735 double a7 = 0.125006;
736 /* JJV changed the initial values of MUPREV and MUOLD */
737 double fact[10] = {
738 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0
739 };
740 /* JJV added ll to the list, for Case A */
741 int32_t ignpoi_mt, j, k, kflag, l, ll, m;
742 double b1, b2, c, c0, c1, c2, c3, d, del, difmuk, e, fk, fx, fy, g, omega, p, p0, px, py, q, s, t, u, v, x, xx, pp[35];
743
744 if (mu < 10.0) goto S120;
745 /*
746 C A S E A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED)
747 JJV changed l in Case A to ll
748 */
749 s = sqrt(mu);
750 d = 6.0 * mu * mu;
751 /*
752 THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
753 PROBABILITIES FK WHENEVER K >= M(MU). LL=IFIX(MU-1.1484)
754 IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
755 */
756 ll = (int32_t)(mu - 1.1484);
757
758 /*
759 STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
760 */
761 g = mu + s * snorm_mt(tn);
762 if (g < 0.0) goto S20;
763 ignpoi_mt = (int32_t)(g);
764 /*
765 STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
766 */
767 if (ignpoi_mt >= ll) return ignpoi_mt;
768 /*
769 STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
770 */
771 fk = (double)ignpoi_mt;
772 difmuk = mu - fk;
773 u = ranf_mt(tn);
774 if (d * u >= difmuk * difmuk * difmuk) return ignpoi_mt;
775S20:
776 /*
777 STEP P. PREPARATIONS FOR STEPS Q AND H.
778 (RECALCULATIONS OF PARAMETERS IF NECESSARY)
779 .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7.
780 THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
781 APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
782 C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
783 */
784 omega = 0.3989423 / s;
785 b1 = 4.166667E-2 / mu;
786 b2 = 0.3 * b1 * b1;
787 c3 = 0.1428571 * b1 * b2;
788 c2 = b2 - 15.0 * c3;
789 c1 = b1 - 6.0 * b2 + 45.0 * c3;
790 c0 = 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
791 c = 0.1069 / mu;
792
793 if (g < 0.0) goto S50;
794 /*
795 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
796 */
797 kflag = 0;
798 goto S70;
799S40:
800 /*
801 STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
802 */
803 if (fy - u * fy <= py * exp(px - fx)) return ignpoi_mt;
804S50:
805 /*
806 STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
807 DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
808 (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
809 */
810 e = sexpo_mt(tn);
811 u = ranf_mt(tn);
812 u += (u - 1.0);
813 t = 1.8 + fsign(e, u);
814 if (t <= -0.6744) goto S50;
815 ignpoi_mt = (int32_t)(mu + s * t);
816 fk = (double)ignpoi_mt;
817 difmuk = mu - fk;
818 /*
819 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
820 */
821 kflag = 1;
822 goto S70;
823S60:
824 /*
825 STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
826 */
827 if (c * fabs(u) > py * exp(px + e) - fy * exp(fx + e)) goto S50;
828 return ignpoi_mt;
829S70:
830 /*
831 STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
832 CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
833 */
834 if (ignpoi_mt >= 10) goto S80;
835 px = -mu;
836 py = pow(mu, (double)ignpoi_mt) / *(fact + ignpoi_mt);
837 goto S110;
838S80:
839 /*
840 CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
841 A0-A7 FOR ACCURACY WHEN ADVISABLE
842 .8333333E-1=1./12. .3989423=(2*PI)**(-.5)
843 */
844 del = 8.333333E-2 / fk;
845 del -= (4.8 * del * del * del);
846 v = difmuk / fk;
847 if (fabs(v) <= 0.25) goto S90;
848 px = fk * log(1.0 + v) - difmuk - del;
849 goto S100;
850S90:
851 px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v + a0) - del;
852S100:
853 py = 0.3989423 / sqrt(fk);
854S110:
855 x = (0.5 - difmuk) / s;
856 xx = x * x;
857 fx = -0.5 * xx;
858 fy = omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
859 if (kflag <= 0) goto S40;
860 goto S60;
861S120:
862 /*
863 C A S E B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
864 JJV changed MUPREV assignment to initial value
865 */
866 m = std::max(INT32_C(1), (int32_t)(mu));
867 l = 0;
868 p = exp(-mu);
869 q = p0 = p;
870S130:
871 /*
872 STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
873 */
874 u = ranf_mt(tn);
875 ignpoi_mt = 0;
876 if (u <= p0) return ignpoi_mt;
877 /*
878 STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
879 PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
880 (0.458=PP(9) FOR MU=10)
881 */
882 if (l == 0) goto S150;
883 j = 1;
884 if (u > 0.458) j = std::min(l, m);
885 for (k = j; k <= l; k++) {
886 if (u <= *(pp + k - 1)) goto S180;
887 }
888 if (l == 35) goto S130;
889S150:
890 /*
891 STEP C. CREATION OF NEW POISSON PROBABILITIES P
892 AND THEIR CUMULATIVES Q=PP(K)
893 */
894 l += 1;
895 for (k = l; k <= 35; k++) {
896 p = p * mu / (double)k;
897 q += p;
898 *(pp + k - 1) = q;
899 if (u <= q) goto S170;
900 }
901 l = 35;
902 goto S130;
903S170:
904 l = k;
905S180:
906 ignpoi_mt = k;
907 return ignpoi_mt;
908}
909int32_t ignbin_mt(int32_t n, double pp, int tn)
910{
911 /*
912**********************************************************************
913int32_t ignbin_mt(int32_t n,double pp)
914GENerate BINomial random deviate
915Function
916Generates a single random deviate from a binomial
917distribution whose number of trials is N and whose
918probability of an event in each trial is P.
919Arguments
920n --> The number of trials in the binomial distribution
921from which a random deviate is to be generated.
922JJV (N >= 0)
923pp --> The probability of an event in each trial of the
924binomial distribution from which a random deviate
925is to be generated.
926JJV (0.0 <= PP <= 1.0)
927ignbin <-- A random deviate yielding the number of events
928from N independent trials, each of which has
929a probability of event P.
930Method
931This is algorithm BTPE from:
932Kachitvichyanukul, V. and Schmeiser, B. W.
933Binomial Random Variate Generation.
934Communications of the ACM, 31, 2
935(February, 1988) 216.
936**********************************************************************
937SUBROUTINE BTPEC(N,PP,ISEED,JX)
938BINOMIAL RANDOM VARIATE GENERATOR
939MEAN .LT. 30 -- INVERSE CDF
940MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA
941FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE
942(SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
943THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
944BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
945BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE
946RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
947USABLE ALGORITHM.
948REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
949"BINOMIAL RANDOM VARIATE GENERATION,"
950COMMUNICATIONS OF THE ACM, FORTHCOMING
951WRITTEN: SEPTEMBER 1980.
952LAST REVISED: MAY 1985, JULY 1987
953REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER
954GENERATOR
955ARGUMENTS
956N : NUMBER OF BERNOULLI TRIALS (INPUT)
957PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
958ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
959JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)
960VARIABLES
961PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
962NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
963XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
964P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
965FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
966M: INTEGER VALUE OF THE CURRENT MODE
967FM: FLOATING POINT VALUE OF THE CURRENT MODE
968XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
969P1: AREA OF THE TRIANGLE
970C: HEIGHT OF THE PARALLELOGRAMS
971XM: CENTER OF THE TRIANGLE
972XL: LEFT END OF THE TRIANGLE
973XR: RIGHT END OF THE TRIANGLE
974AL: TEMPORARY VARIABLE
975XLL: RATE FOR THE LEFT EXPONENTIAL TAIL
976XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL
977P2: AREA OF THE PARALLELOGRAMS
978P3: AREA OF THE LEFT EXPONENTIAL TAIL
979P4: AREA OF THE RIGHT EXPONENTIAL TAIL
980U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
981FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
982FROM THE REGION
983V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
984(REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
985REJECT THE CANDIDATE VALUE
986IX: INTEGER CANDIDATE VALUE
987X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
988AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
989K: ABSOLUTE VALUE OF (IX-M)
990F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
991ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
992ALSO USED IN THE INVERSE TRANSFORMATION
993R: THE RATIO P/Q
994G: CONSTANT USED IN CALCULATION OF PROBABILITY
995MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
996OF F WHEN IX IS GREATER THAN M
997IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
998CALCULATION OF F WHEN IX IS LESS THAN M
999I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
1000AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
1001YNORM: LOGARITHM OF NORMAL BOUND
1002ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
1003X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
1004USED IN THE FINAL ACCEPT/REJECT TEST
1005QN: PROBABILITY OF NO SUCCESS IN N TRIALS
1006REMARK
1007IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
1008SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME
1009COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
1010ARE NOT INVOLVED.
1011ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
1012GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
1013TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
1014**********************************************************************
1015*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
1016*/
1017/* JJV changed initial values to ridiculous values */
1018 double psave = -1.0E37;
1019 int32_t nsave = -214748365;
1020 int32_t ignbin_mt, i, ix, ix1, k, m, mp, T1;
1021 double al, alv, amaxp, c, f, f1, f2, ffm, fm, g, p, p1, p2, p3, p4, q, qn, r, u, v, w, w2, x, x1,
1022 x2, xl, xll, xlr, xm, xnp, xnpq, xr, ynorm, z, z2;
1023
1024 /*
1025 *****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
1026 JJV added checks to ensure 0.0 <= PP <= 1.0
1027 */
1028 if (pp < 0.0) ERR_CRITICAL("PP < 0.0 in IGNBIN");
1029 if (pp > 1.0) ERR_CRITICAL("PP > 1.0 in IGNBIN");
1030 psave = pp;
1031 p = std::min(psave, 1.0 - psave);
1032 q = 1.0 - p;
1033
1034 /*
1035 JJV added check to ensure N >= 0
1036 */
1037 if (n < 0) ERR_CRITICAL("N < 0 in IGNBIN");
1038 xnp = n * p;
1039 nsave = n;
1040 if (xnp < 30.0) goto S140;
1041 ffm = xnp + p;
1042 m = (int32_t)ffm;
1043 fm = m;
1044 xnpq = xnp * q;
1045 p1 = (int32_t)(2.195 * sqrt(xnpq) - 4.6 * q) + 0.5;
1046 xm = fm + 0.5;
1047 xl = xm - p1;
1048 xr = xm + p1;
1049 c = 0.134 + 20.5 / (15.3 + fm);
1050 al = (ffm - xl) / (ffm - xl * p);
1051 xll = al * (1.0 + 0.5 * al);
1052 al = (xr - ffm) / (xr * q);
1053 xlr = al * (1.0 + 0.5 * al);
1054 p2 = p1 * (1.0 + c + c);
1055 p3 = p2 + c / xll;
1056 p4 = p3 + c / xlr;
1057S30:
1058 /*
1059 *****GENERATE VARIATE
1060 */
1061 u = ranf_mt(tn) * p4;
1062 v = ranf_mt(tn);
1063 /*
1064 TRIANGULAR REGION
1065 */
1066 if (u > p1) goto S40;
1067 ix = (int32_t)(xm - p1 * v + u);
1068 goto S170;
1069S40:
1070 /*
1071 PARALLELOGRAM REGION
1072 */
1073 if (u > p2) goto S50;
1074 x = xl + (u - p1) / c;
1075 v = v * c + 1.0 - std::abs(xm - x) / p1;
1076 if (v > 1.0 || v <= 0.0) goto S30;
1077 ix = (int32_t)x;
1078 goto S70;
1079S50:
1080 /*
1081 LEFT TAIL
1082 */
1083 if (u > p3) goto S60;
1084 ix = (int32_t)(xl + log(v) / xll);
1085 if (ix < 0) goto S30;
1086 v *= ((u - p2) * xll);
1087 goto S70;
1088S60:
1089 /*
1090 RIGHT TAIL
1091 */
1092 ix = (int32_t)(xr - log(v) / xlr);
1093 if (ix > n) goto S30;
1094 v *= ((u - p3) * xlr);
1095S70:
1096 /*
1097 *****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
1098 */
1099 k = std::abs(ix - m);
1100 if (k > 20 && k < xnpq / 2 - 1) goto S130;
1101 /*
1102 EXPLICIT EVALUATION
1103 */
1104 f = 1.0;
1105 r = p / q;
1106 g = (n + 1) * r;
1107 T1 = m - ix;
1108 if (T1 < 0) goto S80;
1109 else if (T1 == 0) goto S120;
1110 else goto S100;
1111S80:
1112 mp = m + 1;
1113 for (i = mp; i <= ix; i++) f *= (g / i - r);
1114 goto S120;
1115S100:
1116 ix1 = ix + 1;
1117 for (i = ix1; i <= m; i++) f /= (g / i - r);
1118S120:
1119 if (v <= f) goto S170;
1120 goto S30;
1121S130:
1122 /*
1123 SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
1124 */
1125 amaxp = k / xnpq * ((k * (k / 3.0 + 0.625) + 0.1666666666666) / xnpq + 0.5);
1126 ynorm = -(k * k / (2.0 * xnpq));
1127 alv = log(v);
1128 if (alv < ynorm - amaxp) goto S170;
1129 if (alv > ynorm + amaxp) goto S30;
1130 /*
1131 STIRLING'S FORMULA TO MACHINE ACCURACY FOR
1132 THE FINAL ACCEPTANCE/REJECTION TEST
1133 */
1134 x1 = ix + 1.0;
1135 f1 = fm + 1.0;
1136 z = n + 1.0 - fm;
1137 w = n - ix + 1.0;
1138 z2 = z * z;
1139 x2 = x1 * x1;
1140 f2 = f1 * f1;
1141 w2 = w * w;
1142 if (alv <= xm * log(f1 / x1) + (n - m + 0.5) * log(z / w) + (ix - m) * log(w * p / (x1 * q)) + (13860.0 -
1143 (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 / 166320.0 + (13860.0 - (462.0 -
1144 (132.0 - (99.0 - 140.0 / z2) / z2) / z2) / z2) / z / 166320.0 + (13860.0 - (462.0 - (132.0 -
1145 (99.0 - 140.0 / x2) / x2) / x2) / x2) / x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0
1146 - 140.0 / w2) / w2) / w2) / w2) / w / 166320.0) goto S170;
1147 goto S30;
1148S140:
1149 /*
1150 INVERSE CDF LOGIC FOR MEAN LESS THAN 30
1151 */
1152 qn = pow(q, (double)n);
1153 r = p / q;
1154 g = r * (n + 1);
1155S150:
1156 ix = 0;
1157 f = qn;
1158 u = ranf_mt(tn);
1159S160:
1160 if (u < f) goto S170;
1161 if (ix > 110) goto S150;
1162 u -= f;
1163 ix += 1;
1164 f *= (g / ix - r);
1165 goto S160;
1166S170:
1167 if (psave > 0.5) ix = n - ix;
1168 ignbin_mt = ix;
1169 return ignbin_mt;
1170}
1171double sexpo(void)
1172/*
1173**********************************************************************
1174
1175
1176 (STANDARD-) E X P O N E N T I A L DISTRIBUTION
1177
1178
1179**********************************************************************
1180**********************************************************************
1181
1182 FOR DETAILS SEE:
1183
1184 AHRENS, J.H. AND DIETER, U.
1185 COMPUTER METHODS FOR SAMPLING FROM THE
1186 EXPONENTIAL AND NORMAL DISTRIBUTIONS.
1187 COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
1188
1189 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
1190 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1191
1192 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1193 SUNIF. The argument IR thus goes away.
1194
1195**********************************************************************
1196 Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
1197 (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
1198*/
1199{
1200 static double q[8] = {
1201 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,
1202 .9999999
1203 };
1204 int32_t i;
1205 double sexpo, a, u, ustar, umin;
1206
1207 a = 0.0;
1208 u = ranf();
1209 goto S30;
1210S20:
1211 a += q[0];
1212S30:
1213 u += u;
1214 /*
1215 * JJV changed the following to reflect the true algorithm and prevent
1216 * JJV unpredictable behavior if U is initially 0.5.
1217 * if(u <= 1.0) goto S20;
1218 */
1219 if (u < 1.0) goto S20;
1220 u -= 1.0;
1221 if (u > q[0]) goto S60;
1222 sexpo = a + u;
1223 return sexpo;
1224S60:
1225 i = 1;
1226 ustar = ranf();
1227 umin = ustar;
1228S70:
1229 ustar = ranf();
1230 if (ustar < umin) umin = ustar;
1231 i += 1;
1232 if (u > q[i - 1]) goto S70;
1233 return a + umin * q[0];
1234}
1235
1236double sexpo_mt(int tn)
1237/*
1238**********************************************************************
1239
1240
1241(STANDARD-) E X P O N E N T I A L DISTRIBUTION
1242
1243
1244**********************************************************************
1245**********************************************************************
1246
1247FOR DETAILS SEE:
1248
1249AHRENS, J.H. AND DIETER, U.
1250COMPUTER METHODS FOR SAMPLING FROM THE
1251EXPONENTIAL AND NORMAL DISTRIBUTIONS.
1252COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
1253
1254ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
1255'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1256
1257Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1258SUNIF. The argument IR thus goes away.
1259
1260**********************************************************************
1261Q(N) = SUM(ALOG(2.0)**K/K!) K=1,..,N , THE HIGHEST N
1262(HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
1263*/
1264{
1265// return -log(1 - ranf_mt(tn)); // a much simpler exponential generator!
1266
1267 double q[8] = {0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,0.99999999999999989};
1268 int32_t i;
1269 double sexpo_mt, a, u, ustar, umin;
1270
1271 a = 0.0;
1272 u = ranf_mt(tn);
1273 goto S30;
1274S20:
1275 a += q[0];
1276S30:
1277 u += u;
1278
1279 if (u < 1.0) goto S20;
1280 u -= 1.0;
1281 if (u > q[0]) goto S60;
1282 sexpo_mt = a + u;
1283 return sexpo_mt;
1284S60:
1285 i = 1;
1286 ustar = ranf_mt(tn);
1287 umin = ustar;
1288S70:
1289 ustar = ranf_mt(tn);
1290 if (ustar < umin) umin = ustar;
1291 i += 1;
1292 if (u > q[i - 1]) goto S70;
1293 return a + umin * q[0];
1294
1295}
1296
1297double snorm(void)
1298/*
1299**********************************************************************
1300
1301
1302 (STANDARD-) N O R M A L DISTRIBUTION
1303
1304
1305**********************************************************************
1306**********************************************************************
1307
1308 FOR DETAILS SEE:
1309
1310 AHRENS, J.H. AND DIETER, U.
1311 EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
1312 SAMPLING FROM THE NORMAL DISTRIBUTION.
1313 MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
1314
1315 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
1316 (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1317
1318 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1319 SUNIF. The argument IR thus goes away.
1320
1321**********************************************************************
1322 THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
1323 H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
1324*/
1325{
1326 static double a[32] = {
1327 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
1328 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
1329 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
1330 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
1331 1.862732,2.153875
1332 };
1333 static double d[31] = {
1334 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
1335 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
1336 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
1337 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
1338 };
1339 static double t[31] = {
1340 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
1341 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
1342 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
1343 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
1344 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
1345 };
1346 static double h[31] = {
1347 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
1348 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
1349 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
1350 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
1351 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
1352 };
1353 int32_t i; //made this non-static: ggilani 27/11/14
1354 double snorm, u, s, ustar, aa, w, y, tt; //made this non-static: ggilani 27/11/14
1355 u = ranf();
1356 s = 0.0;
1357 if (u > 0.5) s = 1.0;
1358 u += (u - s);
1359 u = 32.0 * u;
1360 i = (int32_t)(u);
1361 if (i == 32) i = 31;
1362 if (i == 0) goto S100;
1363 /*
1364 START CENTER
1365 */
1366 ustar = u - (double)i;
1367 aa = *(a + i - 1);
1368S40:
1369 if (ustar <= *(t + i - 1)) goto S60;
1370 w = (ustar - *(t + i - 1)) * *(h + i - 1);
1371S50:
1372 /*
1373 EXIT (BOTH CASES)
1374 */
1375 y = aa + w;
1376 snorm = y;
1377 if (s == 1.0) snorm = -y;
1378 return snorm;
1379S60:
1380 /*
1381 CENTER CONTINUED
1382 */
1383 u = ranf();
1384 w = u * (*(a + i) - aa);
1385 tt = (0.5 * w + aa) * w;
1386 goto S80;
1387S70:
1388 tt = u;
1389 ustar = ranf();
1390S80:
1391 if (ustar > tt) goto S50;
1392 u = ranf();
1393 if (ustar >= u) goto S70;
1394 ustar = ranf();
1395 goto S40;
1396S100:
1397 /*
1398 START TAIL
1399 */
1400 i = 6;
1401 aa = *(a + 31);
1402 goto S120;
1403S110:
1404 aa += *(d + i - 1);
1405 i += 1;
1406S120:
1407 u += u;
1408 if (u < 1.0) goto S110;
1409 u -= 1.0;
1410S140:
1411 w = u * *(d + i - 1);
1412 tt = (0.5 * w + aa) * w;
1413 goto S160;
1414S150:
1415 tt = u;
1416S160:
1417 ustar = ranf();
1418 if (ustar > tt) goto S50;
1419 u = ranf();
1420 if (ustar >= u) goto S150;
1421 u = ranf();
1422 goto S140;
1423}
1424double snorm_mt(int tn)
1425/*
1426**********************************************************************
1427
1428
1429(STANDARD-) N O R M A L DISTRIBUTION
1430
1431
1432**********************************************************************
1433**********************************************************************
1434
1435FOR DETAILS SEE:
1436
1437AHRENS, J.H. AND DIETER, U.
1438EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
1439SAMPLING FROM THE NORMAL DISTRIBUTION.
1440MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
1441
1442ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
1443(M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1444
1445Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1446SUNIF. The argument IR thus goes away.
1447
1448**********************************************************************
1449THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
1450H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
1451*/
1452{
1453 static double a[32] = {
1454 0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
1455 0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
1456 0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
1457 1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
1458 1.862732,2.153875
1459 };
1460 static double d[31] = {
1461 0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
1462 0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
1463 0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
1464 0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
1465 };
1466 static double t[31] = {
1467 7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
1468 1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
1469 2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
1470 4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
1471 9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
1472 };
1473 static double h[31] = {
1474 3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
1475 4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
1476 4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
1477 5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
1478 8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
1479 };
1480 int32_t i;
1481 double snorm_mt, u, s, ustar, aa, w, y, tt;
1482 u = ranf_mt(tn);
1483 s = 0.0;
1484 if (u > 0.5) s = 1.0;
1485 u += (u - s);
1486 u = 32.0 * u;
1487 i = (int32_t)(u);
1488 if (i == 32) i = 31;
1489 if (i == 0) goto S100;
1490 /*
1491 START CENTER
1492 */
1493 ustar = u - (double)i;
1494 aa = *(a + i - 1);
1495S40:
1496 if (ustar <= *(t + i - 1)) goto S60;
1497 w = (ustar - *(t + i - 1)) * *(h + i - 1);
1498S50:
1499 /*
1500 EXIT (BOTH CASES)
1501 */
1502 y = aa + w;
1503 snorm_mt = y;
1504 if (s == 1.0) snorm_mt = -y;
1505 return snorm_mt;
1506S60:
1507 /*
1508 CENTER CONTINUED
1509 */
1510 u = ranf_mt(tn);
1511 w = u * (*(a + i) - aa);
1512 tt = (0.5 * w + aa) * w;
1513 goto S80;
1514S70:
1515 tt = u;
1516 ustar = ranf_mt(tn);
1517S80:
1518 if (ustar > tt) goto S50;
1519 u = ranf_mt(tn);
1520 if (ustar >= u) goto S70;
1521 ustar = ranf_mt(tn);
1522 goto S40;
1523S100:
1524 /*
1525 START TAIL
1526 */
1527 i = 6;
1528 aa = *(a + 31);
1529 goto S120;
1530S110:
1531 aa += *(d + i - 1);
1532 i += 1;
1533S120:
1534 u += u;
1535 if (u < 1.0) goto S110;
1536 u -= 1.0;
1537S140:
1538 w = u * *(d + i - 1);
1539 tt = (0.5 * w + aa) * w;
1540 goto S160;
1541S150:
1542 tt = u;
1543S160:
1544 ustar = ranf_mt(tn);
1545 if (ustar > tt) goto S50;
1546 u = ranf_mt(tn);
1547 if (ustar >= u) goto S150;
1548 u = ranf_mt(tn);
1549 goto S140;
1550}
1551double fsign(double num, double sign)
1552/* Transfers sign of argument sign to argument num */
1553{
1554 if ((sign > 0.0f && num < 0.0f) || (sign < 0.0f && num>0.0f))
1555 return -num;
1556 else return num;
1557}
1558/*function gen_snorm
1559 * purpose: my own implementation of sampling from a uniform distribution, using Marsaglia polar method, but for multi-threading
1560 *
1561 * author: ggilani, date: 28/11/14
1562 */
1563double gen_norm_mt(double mu, double sd, int tn)
1564{
1565 double u, v, x, S;
1566
1567 do
1568 {
1569 u = 2 * ranf_mt(tn) - 1; //u and v are uniform random numbers on the interval [-1,1]
1570 v = 2 * ranf_mt(tn) - 1;
1571
1572 //calculate S=U^2+V^2
1573 S = u * u + v * v;
1574 } while (S >= 1 || S == 0);
1575
1576 //calculate x,y - both of which are normally distributed
1577 x = u * sqrt((-2 * log(S)) / S);
1578
1579 // This routine can be accelerated by storing the second normal value
1580 // and using it for the next call.
1581 // y = v * sqrt((-2 * log(S)) / S);
1582
1583 //return x
1584 return x * sd + mu;
1585}
1586
1587/*function gen_gamma_mt
1588 * purpose: my own implementation of sampling from a gamma distribution, using Marsaglia-Tsang method, but for multi-threading
1589 *
1590 * author: ggilani, date: 01/12/14
1591 */
1592double gen_gamma_mt(double beta, double alpha, int tn)
1593{
1594 double d, c, u, v, z, f, alpha2, gamma;
1595
1596 //error statment if either beta or alpha are <=0, as gamma distribution is undefined in this case
1597 if ((beta <= 0) || (alpha <= 0))
1598 {
1599 ERR_CRITICAL("Gamma distribution parameters in undefined range!\n");
1600 }
1601
1602 //method is slightly different depending on whether alpha is greater than or equal to 1, or less than 1
1603 if (alpha >= 1)
1604 {
1605 d = alpha - (1.0 / 3.0);
1606 c = 1.0 / (3.0 * sqrt(d));
1607 do
1608 {
1609 //sample one random number from uniform distribution and one from standard normal distribution
1610 u = ranf_mt(tn);
1611 z = gen_norm_mt(0, 1, tn);
1612 v = 1 + z * c;
1613 v = v * v * v;
1614 } while ((z <= (-1.0 / c)) || (log(u) >= (0.5 * z * z + d - d * v + d * log(v))));
1615 //if beta is equal to 1, there is no scale. If beta is not equal to one, return scaled value
1616 if (beta != 1)
1617 {
1618 return (d * v) / beta;
1619 }
1620 else
1621 {
1622 return d * v;
1623 }
1624 }
1625 else
1626 {
1627 //if alpha is less than 1, initially sample from gamma(beta,alpha+1)
1628 alpha2 = alpha + 1;
1629 d = alpha2 - (1.0 / 3.0);
1630 c = 1.0 / (3.0 * sqrt(d));
1631 do
1632 {
1633 //sample one random number from uniform distribution and one from standard normal distribution
1634 u = ranf_mt(tn);
1635 z = gen_norm_mt(0, 1, tn);
1636 v = 1 + z * c;
1637 v = v * v * v;
1638 } while ((z <= (-1.0 / c)) || (log(u) >= (0.5 * z * z + d - d * v + d * log(v))));
1639 //if beta is equal to 1, there is no scale. If beta is not equal to one, return scaled value
1640 if (beta != 1)
1641 {
1642 gamma = (d * v) / beta;
1643 }
1644 else
1645 {
1646 gamma = d * v;
1647 }
1648 //now rescale again to take into account that alpha is less than 1
1649 f = pow(ranf_mt(tn), (1.0 / alpha));
1650 //return gamma scaled by f
1651 return gamma * f;
1652 }
1653}
1654/* function gen_lognormal(double mu, double sigma)
1655 * purpose: to generate samples from a lognormal distribution with parameters mu and sigma
1656 *
1657 * parameters:
1658 * mean mu
1659 * standard deviation sigma
1660 *
1661 * returns: double from the specified lognormal distribution
1662 *
1663 * author: ggilani, date: 09/02/17
1664 */
1665double gen_lognormal(double mu, double sigma)
1666{
1667 double randnorm, location, scale;
1668
1669 randnorm = snorm();
1670 location = log(mu / sqrt(1 + ((sigma * sigma) / (mu * mu))));
1671 scale = sqrt(log(1 + ((sigma * sigma) / (mu * mu))));
1672 return exp(location + scale * randnorm);
1673}
1674void SampleWithoutReplacement(int tn, int k, int n)
1675{
1676 /* Based on algorithm SG of http://portal.acm.org/citation.cfm?id=214402
1677 ACM Transactions on Mathematical Software (TOMS) archive
1678 Volume 11 , Issue 2 (June 1985) table of contents
1679 Pages: 157 - 169
1680 Year of Publication: 1985
1681 ISSN:0098-3500
1682 */
1683
1684 double t, r, a, mu, f;
1685 int i, j, q, b;
1686
1687 if (k < 3)
1688 {
1689 for (i = 0; i < k; i++)
1690 {
1691 do
1692 {
1693 SamplingQueue[tn][i] = (int)(ranf_mt(tn) * ((double)n));
1694// This original formulation is completely valid, but the PVS Studio analyzer
1695// notes this, so I am changing it just to get report-clean.
1696// "V1008 Consider inspecting the 'for' operator. No more than one iteration of the loop will be performed. Rand.cpp 2450"
1697// for (j = q = 0; (j < i) && (!q); j++)
1698// q = (SamplingQueue[tn][i] == SamplingQueue[tn][j]);
1699 j = q = 0;
1700 if (i == 1)
1701 q = (SamplingQueue[tn][i] == SamplingQueue[tn][j]);
1702 } while (q);
1703 }
1704 q = k;
1705 }
1706 else if (2 * k > n)
1707 {
1708 for (i = 0; i < n; i++)
1709 SamplingQueue[tn][i] = i;
1710 for (i = n; i > k; i--)
1711 {
1712 j = (int)(ranf_mt(tn) * ((double)i));
1713 if (j != i - 1)
1714 {
1715 b = SamplingQueue[tn][j];
1716 SamplingQueue[tn][j] = SamplingQueue[tn][i - 1];
1717 SamplingQueue[tn][i - 1] = b;
1718 }
1719 }
1720 q = k;
1721 }
1722 else if (4 * k > n)
1723 {
1724 for (i = 0; i < n; i++)
1725 SamplingQueue[tn][i] = i;
1726 for (i = 0; i < k; i++)
1727 {
1728 j = (int)(ranf_mt(tn) * ((double)(n - i)));
1729 if (j > 0)
1730 {
1731 b = SamplingQueue[tn][i];
1732 SamplingQueue[tn][i] = SamplingQueue[tn][i + j];
1733 SamplingQueue[tn][i + j] = b;
1734 }
1735 }
1736 q = k;
1737 }
1738 else
1739 {
1740 /* fprintf(stderr,"@%i %i:",k,n); */
1741 t = (double)k;
1742 r = sqrt(t);
1743 a = sqrt(log(1 + t / 2 * PI));
1744 a = a + a * a / (3 * r);
1745 mu = t + a * r;
1746 b = 2 * MAX_PLACE_SIZE; /* (int) (k+4*a*r); */
1747 f = -1 / (log(1 - mu / ((double)n)));
1748 i = q = 0;
1749 while (i <= n)
1750 {
1751 i += (int)ceil(-log(ranf_mt(tn)) * f);
1752 if (i <= n)
1753 {
1754 SamplingQueue[tn][q] = i - 1;
1755 q++;
1756 if (q >= b) i = q = 0;
1757 }
1758 else if (q < k)
1759 i = q = 0;
1760 }
1761 }
1762 /* else
1763 {
1764 t=(double) (n-k);
1765 r=sqrt(t);
1766 a=sqrt(log(1+t/2*PI));
1767 a=a+a*a/(3*r);
1768 mu=t+a*r;
1769 b=2*MAX_PLACE_SIZE;
1770 f=-1/(log(1-mu/((double) n)));
1771 i=q=0;
1772 while(i<=n)
1773 {
1774 int i2=i+(int) ceil(-log(ranf_mt(tn))*f);
1775 i++;
1776 if(i2<=n)
1777 for(;(i<i2)&&(q<b);i++)
1778 {
1779 SamplingQueue[tn][q]=i-1;
1780 q++;
1781 }
1782 else
1783 {
1784 for(;(i<=n)&&(q<b);i++)
1785 {
1786 SamplingQueue[tn][q]=i-1;
1787 q++;
1788 }
1789 if(q<k) i=q=0;
1790 }
1791 if(q>=b) i=q=0;
1792 }
1793 }
1794 */
1795 /* if(k>2)
1796 {
1797 fprintf(stderr,"(%i) ",q);
1798 for(i=0;i<q;i++) fprintf(stderr,"%i ",SamplingQueue[tn][i]);
1799 fprintf(stderr,"\n");
1800 }
1801 */ while (q > k)
1802 {
1803 i = (int)(ranf_mt(tn) * ((double)q));
1804 if (i < q - 1) SamplingQueue[tn][i] = SamplingQueue[tn][q - 1];
1805 q--;
1806 }
1807
1808}
1809
Generated by doxygen 1.14.0