covid-sim
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 */
13 int32_t* Xcg1, *Xcg2;
14 int **SamplingQueue = nullptr;
15 
19 
20 double 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 
45 double 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 
66 void 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 
100 int32_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 
177 int32_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;
325 S30:
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;
337 S40:
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;
347 S50:
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;
356 S60:
357  /*
358  RIGHT TAIL
359  */
360  ix = (int32_t)(xr - log(v) / xlr);
361  if (ix > n) goto S30;
362  v *= ((u - p3) * xlr);
363 S70:
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;
379 S80:
380  mp = m + 1;
381  for (i = mp; i <= ix; i++) f *= (g / i - r);
382  goto S120;
383 S100:
384  ix1 = ix + 1;
385  for (i = ix1; i <= m; i++) f /= (g / i - r);
386 S120:
387  if (v <= f) goto S170;
388  goto S30;
389 S130:
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;
416 S140:
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);
423 S150:
424  ix = 0;
425  f = qn;
426  u = ranf();
427 S160:
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;
434 S170:
435  if (psave > 0.5) ix = n - ix;
436  ignbin = ix;
437  return ignbin;
438 }
439 
440 int32_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;
540 S20:
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;
564 S40:
565  /*
566  STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
567  */
568  if (fy - u * fy <= py * exp(px - fx)) return ignpoi;
569 S50:
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;
588 S60:
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;
594 S70:
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;
603 S80:
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;
615 S90:
616  px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v + a0) - del;
617 S100:
618  py = 0.3989423 / sqrt(fk);
619 S110:
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;
626 S120:
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;
635 S130:
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;
654 S150:
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;
668 S170:
669  l = k;
670 S180:
671  ignpoi = k;
672  return ignpoi;
673 }
674 
675 int32_t ignpoi_mt(double mu, int tn)
676 /*
677 **********************************************************************
678 int32_t ignpoi_mt(double mu)
679 GENerate POIsson random deviate
680 Function
681 Generates a single random deviate from a Poisson
682 distribution with mean MU.
683 Arguments
684 mu --> The mean of the Poisson distribution from which
685 a random deviate is to be generated.
686 (mu >= 0.0)
687 ignpoi_mt <-- The random deviate.
688 Method
689 Renames KPOIS from TOMS as slightly modified by BWB to use RANF
690 instead of SUNIF.
691 For details see:
692 Ahrens, J.H. and Dieter, U.
693 Computer Generation of Poisson Deviates
694 From Modified Normal Distributions.
695 ACM Trans. Math. Software, 8, 2
696 (June 1982),163-179
697 **********************************************************************
698 **********************************************************************
699 
700 
701 P O I S S O N DISTRIBUTION
702 
703 
704 **********************************************************************
705 **********************************************************************
706 
707 FOR DETAILS SEE:
708 
709 AHRENS, J.H. AND DIETER, U.
710 COMPUTER GENERATION OF POISSON DEVIATES
711 FROM MODIFIED NORMAL DISTRIBUTIONS.
712 ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
713 
714 (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
715 
716 **********************************************************************
717 INTEGER FUNCTION IGNPOI(IR,MU)
718 INPUT: IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
719 MU=MEAN MU OF THE POISSON DISTRIBUTION
720 OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
721 MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
722 TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
723 COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
724 SEPARATION 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;
775 S20:
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;
799 S40:
800  /*
801  STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
802  */
803  if (fy - u * fy <= py * exp(px - fx)) return ignpoi_mt;
804 S50:
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;
823 S60:
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;
829 S70:
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;
838 S80:
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;
850 S90:
851  px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v + a0) - del;
852 S100:
853  py = 0.3989423 / sqrt(fk);
854 S110:
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;
861 S120:
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;
870 S130:
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;
889 S150:
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;
903 S170:
904  l = k;
905 S180:
906  ignpoi_mt = k;
907  return ignpoi_mt;
908 }
909 int32_t ignbin_mt(int32_t n, double pp, int tn)
910 {
911  /*
912 **********************************************************************
913 int32_t ignbin_mt(int32_t n,double pp)
914 GENerate BINomial random deviate
915 Function
916 Generates a single random deviate from a binomial
917 distribution whose number of trials is N and whose
918 probability of an event in each trial is P.
919 Arguments
920 n --> The number of trials in the binomial distribution
921 from which a random deviate is to be generated.
922 JJV (N >= 0)
923 pp --> The probability of an event in each trial of the
924 binomial distribution from which a random deviate
925 is to be generated.
926 JJV (0.0 <= PP <= 1.0)
927 ignbin <-- A random deviate yielding the number of events
928 from N independent trials, each of which has
929 a probability of event P.
930 Method
931 This is algorithm BTPE from:
932 Kachitvichyanukul, V. and Schmeiser, B. W.
933 Binomial Random Variate Generation.
934 Communications of the ACM, 31, 2
935 (February, 1988) 216.
936 **********************************************************************
937 SUBROUTINE BTPEC(N,PP,ISEED,JX)
938 BINOMIAL RANDOM VARIATE GENERATOR
939 MEAN .LT. 30 -- INVERSE CDF
940 MEAN .GE. 30 -- ALGORITHM BTPE: ACCEPTANCE-REJECTION VIA
941 FOUR REGION COMPOSITION. THE FOUR REGIONS ARE A TRIANGLE
942 (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
943 THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
944 BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
945 BTPEC REFERS TO BTPE AND "COMBINED." THUS BTPE IS THE
946 RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
947 USABLE ALGORITHM.
948 REFERENCE: VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
949 "BINOMIAL RANDOM VARIATE GENERATION,"
950 COMMUNICATIONS OF THE ACM, FORTHCOMING
951 WRITTEN: SEPTEMBER 1980.
952 LAST REVISED: MAY 1985, JULY 1987
953 REQUIRED SUBPROGRAM: RAND() -- A UNIFORM (0,1) RANDOM NUMBER
954 GENERATOR
955 ARGUMENTS
956 N : NUMBER OF BERNOULLI TRIALS (INPUT)
957 PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
958 ISEED: RANDOM NUMBER SEED (INPUT AND OUTPUT)
959 JX: RANDOMLY GENERATED OBSERVATION (OUTPUT)
960 VARIABLES
961 PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
962 NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
963 XNP: VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
964 P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
965 FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
966 M: INTEGER VALUE OF THE CURRENT MODE
967 FM: FLOATING POINT VALUE OF THE CURRENT MODE
968 XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
969 P1: AREA OF THE TRIANGLE
970 C: HEIGHT OF THE PARALLELOGRAMS
971 XM: CENTER OF THE TRIANGLE
972 XL: LEFT END OF THE TRIANGLE
973 XR: RIGHT END OF THE TRIANGLE
974 AL: TEMPORARY VARIABLE
975 XLL: RATE FOR THE LEFT EXPONENTIAL TAIL
976 XLR: RATE FOR THE RIGHT EXPONENTIAL TAIL
977 P2: AREA OF THE PARALLELOGRAMS
978 P3: AREA OF THE LEFT EXPONENTIAL TAIL
979 P4: AREA OF THE RIGHT EXPONENTIAL TAIL
980 U: A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
981 FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
982 FROM THE REGION
983 V: A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
984 (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
985 REJECT THE CANDIDATE VALUE
986 IX: INTEGER CANDIDATE VALUE
987 X: PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
988 AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
989 K: ABSOLUTE VALUE OF (IX-M)
990 F: THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
991 ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
992 ALSO USED IN THE INVERSE TRANSFORMATION
993 R: THE RATIO P/Q
994 G: CONSTANT USED IN CALCULATION OF PROBABILITY
995 MP: MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
996 OF F WHEN IX IS GREATER THAN M
997 IX1: CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
998 CALCULATION OF F WHEN IX IS LESS THAN M
999 I: INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
1000 AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
1001 YNORM: LOGARITHM OF NORMAL BOUND
1002 ALV: NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
1003 X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
1004 USED IN THE FINAL ACCEPT/REJECT TEST
1005 QN: PROBABILITY OF NO SUCCESS IN N TRIALS
1006 REMARK
1007 IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
1008 SAVE A MEMORY POSITION AND A LINE OF CODE. HOWEVER, SOME
1009 COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
1010 ARE NOT INVOLVED.
1011 ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
1012 GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
1013 TYPE 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;
1057 S30:
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;
1069 S40:
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;
1079 S50:
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;
1088 S60:
1089  /*
1090  RIGHT TAIL
1091  */
1092  ix = (int32_t)(xr - log(v) / xlr);
1093  if (ix > n) goto S30;
1094  v *= ((u - p3) * xlr);
1095 S70:
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;
1111 S80:
1112  mp = m + 1;
1113  for (i = mp; i <= ix; i++) f *= (g / i - r);
1114  goto S120;
1115 S100:
1116  ix1 = ix + 1;
1117  for (i = ix1; i <= m; i++) f /= (g / i - r);
1118 S120:
1119  if (v <= f) goto S170;
1120  goto S30;
1121 S130:
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;
1148 S140:
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);
1155 S150:
1156  ix = 0;
1157  f = qn;
1158  u = ranf_mt(tn);
1159 S160:
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;
1166 S170:
1167  if (psave > 0.5) ix = n - ix;
1168  ignbin_mt = ix;
1169  return ignbin_mt;
1170 }
1171 double 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;
1210 S20:
1211  a += q[0];
1212 S30:
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;
1224 S60:
1225  i = 1;
1226  ustar = ranf();
1227  umin = ustar;
1228 S70:
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 
1236 double 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 
1247 FOR DETAILS SEE:
1248 
1249 AHRENS, J.H. AND DIETER, U.
1250 COMPUTER METHODS FOR SAMPLING FROM THE
1251 EXPONENTIAL AND NORMAL DISTRIBUTIONS.
1252 COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
1253 
1254 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
1255 'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1256 
1257 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1258 SUNIF. The argument IR thus goes away.
1259 
1260 **********************************************************************
1261 Q(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;
1274 S20:
1275  a += q[0];
1276 S30:
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;
1284 S60:
1285  i = 1;
1286  ustar = ranf_mt(tn);
1287  umin = ustar;
1288 S70:
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 
1297 double 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);
1368 S40:
1369  if (ustar <= *(t + i - 1)) goto S60;
1370  w = (ustar - *(t + i - 1)) * *(h + i - 1);
1371 S50:
1372  /*
1373  EXIT (BOTH CASES)
1374  */
1375  y = aa + w;
1376  snorm = y;
1377  if (s == 1.0) snorm = -y;
1378  return snorm;
1379 S60:
1380  /*
1381  CENTER CONTINUED
1382  */
1383  u = ranf();
1384  w = u * (*(a + i) - aa);
1385  tt = (0.5 * w + aa) * w;
1386  goto S80;
1387 S70:
1388  tt = u;
1389  ustar = ranf();
1390 S80:
1391  if (ustar > tt) goto S50;
1392  u = ranf();
1393  if (ustar >= u) goto S70;
1394  ustar = ranf();
1395  goto S40;
1396 S100:
1397  /*
1398  START TAIL
1399  */
1400  i = 6;
1401  aa = *(a + 31);
1402  goto S120;
1403 S110:
1404  aa += *(d + i - 1);
1405  i += 1;
1406 S120:
1407  u += u;
1408  if (u < 1.0) goto S110;
1409  u -= 1.0;
1410 S140:
1411  w = u * *(d + i - 1);
1412  tt = (0.5 * w + aa) * w;
1413  goto S160;
1414 S150:
1415  tt = u;
1416 S160:
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 }
1424 double snorm_mt(int tn)
1425 /*
1426 **********************************************************************
1427 
1428 
1429 (STANDARD-) N O R M A L DISTRIBUTION
1430 
1431 
1432 **********************************************************************
1433 **********************************************************************
1434 
1435 FOR DETAILS SEE:
1436 
1437 AHRENS, J.H. AND DIETER, U.
1438 EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM
1439 SAMPLING FROM THE NORMAL DISTRIBUTION.
1440 MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.
1441 
1442 ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'
1443 (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
1444 
1445 Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
1446 SUNIF. The argument IR thus goes away.
1447 
1448 **********************************************************************
1449 THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
1450 H(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);
1495 S40:
1496  if (ustar <= *(t + i - 1)) goto S60;
1497  w = (ustar - *(t + i - 1)) * *(h + i - 1);
1498 S50:
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;
1506 S60:
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;
1514 S70:
1515  tt = u;
1516  ustar = ranf_mt(tn);
1517 S80:
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;
1523 S100:
1524  /*
1525  START TAIL
1526  */
1527  i = 6;
1528  aa = *(a + 31);
1529  goto S120;
1530 S110:
1531  aa += *(d + i - 1);
1532  i += 1;
1533 S120:
1534  u += u;
1535  if (u < 1.0) goto S110;
1536  u -= 1.0;
1537 S140:
1538  w = u * *(d + i - 1);
1539  tt = (0.5 * w + aa) * w;
1540  goto S160;
1541 S150:
1542  tt = u;
1543 S160:
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 }
1551 double 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  */
1563 double 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  */
1592 double 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  */
1665 double 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 }
1674 void 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