5 #include "MachineDefines.h" 14 int **SamplingQueue =
nullptr;
26 curntg = CACHE_LINE_SIZE * omp_get_thread_num();
33 s1 = Xa1 * (s1 - k * 53668) - k * 12211;
34 if (s1 < 0) s1 += Xm1;
36 s2 = Xa2 * (s2 - k * 52774) - k * 3791;
37 if (s2 < 0) s2 += Xm2;
41 if (z < 1) z += (Xm1 - 1);
42 return ((
double)z) / Xm1;
45 double ranf_mt(
int tn)
50 curntg = CACHE_LINE_SIZE * tn;
54 s1 = Xa1 * (s1 - k * 53668) - k * 12211;
55 if (s1 < 0) s1 += Xm1;
57 s2 = Xa2 * (s2 - k * 52774) - k * 3791;
58 if (s2 < 0) s2 += Xm2;
62 if (z < 1) z += (Xm1 - 1);
63 return ((
double)z) / Xm1;
66 void setall(int32_t *pseed1, int32_t *pseed2)
88 int32_t iseed1 = *pseed1;
89 int32_t iseed2 = *pseed2;
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);
100 int32_t mltmod(int32_t a, int32_t s, int32_t m)
114 const int32_t h = 32768;
115 int32_t a0, a1, k, p, q, qh, rh;
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);
138 p = h * (s - k * qh) - k * rh;
149 p -= (k * (m - a1 * q));
151 p += (a1 * (s - k * q));
158 p = h * (p - k * qh) - k * rh;
167 p -= (k * (m - a0 * q));
169 p += (a0 * (s - k * q));
177 int32_t ignbin(int32_t n,
double pp)
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;
296 if (pp < 0.0F) ERR_CRITICAL(
"PP < 0.0 in IGNBIN");
297 if (pp > 1.0F) ERR_CRITICAL(
"PP > 1.0 in IGNBIN");
299 p = std::min(psave, 1.0 - psave);
305 if (n < 0L) ERR_CRITICAL(
"N < 0 in IGNBIN");
308 if (xnp < 30.0)
goto S140;
313 p1 = (int32_t)(2.195 * sqrt(xnpq) - 4.6 * q) + 0.5;
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);
334 if (u > p1)
goto S40;
335 ix = (int32_t)(xm - p1 * v + u);
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;
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);
360 ix = (int32_t)(xr - log(v) / xlr);
361 if (ix > n)
goto S30;
362 v *= ((u - p3) * xlr);
367 k = std::abs(ix - m);
368 if (k > 20 && k < xnpq / 2 - 1)
goto S130;
376 if (T1 < 0)
goto S80;
377 else if (T1 == 0)
goto S120;
381 for (i = mp; i <= ix; i++) f *= (g / i - r);
385 for (i = ix1; i <= m; i++) f /= (g / i - r);
387 if (v <= f)
goto S170;
393 amaxp = k / xnpq * ((k * (k / 3.0 + 0.625) + 0.1666666666666) / xnpq + 0.5);
394 ynorm = -(k * k / (2.0 * xnpq));
396 if (alv < ynorm - amaxp)
goto S170;
397 if (alv > ynorm + amaxp)
goto S30;
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;
420 qn = pow(q, (
double)n);
428 if (u < f)
goto S170;
429 if (ix > 110)
goto S150;
435 if (psave > 0.5) ix = n - ix;
440 int32_t ignpoi(
double mu)
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;
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
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];
509 if (mu < 10.0)
goto S120;
521 ll = (int32_t)(mu - 1.1484);
526 g = mu + s * snorm();
527 if (g < 0.0)
goto S20;
528 ignpoi = (int32_t)(g);
532 if (ignpoi >= ll)
return ignpoi;
539 if (d * u >= difmuk * difmuk * difmuk)
return ignpoi;
549 omega = 0.3989423 / s;
550 b1 = 4.166667E-2 / mu;
552 c3 = 0.1428571 * b1 * b2;
554 c1 = b1 - 6.0 * b2 + 45.0 * c3;
555 c0 = 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
558 if (g < 0.0)
goto S50;
568 if (fy - u * fy <= py * exp(px - fx))
return ignpoi;
578 t = 1.8 + fsign(e, u);
579 if (t <= -0.6744)
goto S50;
580 ignpoi = (int32_t)(mu + s * t);
592 if (c * fabs(u) > py * exp(px + e) - fy * exp(fx + e))
goto S50;
599 if (ignpoi >= 10)
goto S80;
601 py = pow(mu, (
double)ignpoi) / *(fact + ignpoi);
609 del = 8.333333E-2 / fk;
610 del -= (4.8 * del * del * del);
612 if (fabs(v) <= 0.25)
goto S90;
613 px = fk * log(1.0 + v) - difmuk - del;
616 px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v + a0) - del;
618 py = 0.3989423 / sqrt(fk);
620 x = (0.5 - difmuk) / s;
623 fy = omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
624 if (kflag <= 0)
goto S40;
631 m = std::max(INT32_C(1), (int32_t)(mu));
641 if (u <= p0)
return ignpoi;
647 if (l == 0)
goto S150;
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;
653 if (l == 35)
goto S130;
660 for (k = l; k <= 35; k++) {
661 p = p * mu / (double)k;
664 if (u <= q)
goto S170;
675 int32_t ignpoi_mt(
double mu,
int tn)
727 extern double fsign(
double num,
double sign);
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;
738 1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0
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];
744 if (mu < 10.0)
goto S120;
756 ll = (int32_t)(mu - 1.1484);
761 g = mu + s * snorm_mt(tn);
762 if (g < 0.0)
goto S20;
763 ignpoi_mt = (int32_t)(g);
767 if (ignpoi_mt >= ll)
return ignpoi_mt;
771 fk = (double)ignpoi_mt;
774 if (d * u >= difmuk * difmuk * difmuk)
return ignpoi_mt;
784 omega = 0.3989423 / s;
785 b1 = 4.166667E-2 / mu;
787 c3 = 0.1428571 * b1 * b2;
789 c1 = b1 - 6.0 * b2 + 45.0 * c3;
790 c0 = 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
793 if (g < 0.0)
goto S50;
803 if (fy - u * fy <= py * exp(px - fx))
return ignpoi_mt;
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;
827 if (c * fabs(u) > py * exp(px + e) - fy * exp(fx + e))
goto S50;
834 if (ignpoi_mt >= 10)
goto S80;
836 py = pow(mu, (
double)ignpoi_mt) / *(fact + ignpoi_mt);
844 del = 8.333333E-2 / fk;
845 del -= (4.8 * del * del * del);
847 if (fabs(v) <= 0.25)
goto S90;
848 px = fk * log(1.0 + v) - difmuk - del;
851 px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v + a0) - del;
853 py = 0.3989423 / sqrt(fk);
855 x = (0.5 - difmuk) / s;
858 fy = omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
859 if (kflag <= 0)
goto S40;
866 m = std::max(INT32_C(1), (int32_t)(mu));
876 if (u <= p0)
return ignpoi_mt;
882 if (l == 0)
goto S150;
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;
888 if (l == 35)
goto S130;
895 for (k = l; k <= 35; k++) {
896 p = p * mu / (double)k;
899 if (u <= q)
goto S170;
909 int32_t ignbin_mt(int32_t n,
double pp,
int tn)
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;
1028 if (pp < 0.0) ERR_CRITICAL(
"PP < 0.0 in IGNBIN");
1029 if (pp > 1.0) ERR_CRITICAL(
"PP > 1.0 in IGNBIN");
1031 p = std::min(psave, 1.0 - psave);
1037 if (n < 0) ERR_CRITICAL(
"N < 0 in IGNBIN");
1040 if (xnp < 30.0)
goto S140;
1045 p1 = (int32_t)(2.195 * sqrt(xnpq) - 4.6 * q) + 0.5;
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);
1061 u = ranf_mt(tn) * p4;
1066 if (u > p1)
goto S40;
1067 ix = (int32_t)(xm - p1 * v + u);
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;
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);
1092 ix = (int32_t)(xr - log(v) / xlr);
1093 if (ix > n)
goto S30;
1094 v *= ((u - p3) * xlr);
1099 k = std::abs(ix - m);
1100 if (k > 20 && k < xnpq / 2 - 1)
goto S130;
1108 if (T1 < 0)
goto S80;
1109 else if (T1 == 0)
goto S120;
1113 for (i = mp; i <= ix; i++) f *= (g / i - r);
1117 for (i = ix1; i <= m; i++) f /= (g / i - r);
1119 if (v <= f)
goto S170;
1125 amaxp = k / xnpq * ((k * (k / 3.0 + 0.625) + 0.1666666666666) / xnpq + 0.5);
1126 ynorm = -(k * k / (2.0 * xnpq));
1128 if (alv < ynorm - amaxp)
goto S170;
1129 if (alv > ynorm + amaxp)
goto S30;
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;
1152 qn = pow(q, (
double)n);
1160 if (u < f)
goto S170;
1161 if (ix > 110)
goto S150;
1167 if (psave > 0.5) ix = n - ix;
1200 static double q[8] = {
1201 0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,
1205 double sexpo, a, u, ustar, umin;
1219 if (u < 1.0)
goto S20;
1221 if (u > q[0])
goto S60;
1230 if (ustar < umin) umin = ustar;
1232 if (u > q[i - 1])
goto S70;
1233 return a + umin * q[0];
1236 double sexpo_mt(
int tn)
1267 double q[8] = {0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,0.99999999999999989};
1269 double sexpo_mt, a, u, ustar, umin;
1279 if (u < 1.0)
goto S20;
1281 if (u > q[0])
goto S60;
1286 ustar = ranf_mt(tn);
1289 ustar = ranf_mt(tn);
1290 if (ustar < umin) umin = ustar;
1292 if (u > q[i - 1])
goto S70;
1293 return a + umin * q[0];
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,
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
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
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
1354 double snorm, u, s, ustar, aa, w, y, tt;
1357 if (u > 0.5) s = 1.0;
1361 if (i == 32) i = 31;
1362 if (i == 0)
goto S100;
1366 ustar = u - (double)i;
1369 if (ustar <= *(t + i - 1))
goto S60;
1370 w = (ustar - *(t + i - 1)) * *(h + i - 1);
1377 if (s == 1.0) snorm = -y;
1384 w = u * (*(a + i) - aa);
1385 tt = (0.5 * w + aa) * w;
1391 if (ustar > tt)
goto S50;
1393 if (ustar >= u)
goto S70;
1408 if (u < 1.0)
goto S110;
1411 w = u * *(d + i - 1);
1412 tt = (0.5 * w + aa) * w;
1418 if (ustar > tt)
goto S50;
1420 if (ustar >= u)
goto S150;
1424 double snorm_mt(
int tn)
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,
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
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
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
1481 double snorm_mt, u, s, ustar, aa, w, y, tt;
1484 if (u > 0.5) s = 1.0;
1488 if (i == 32) i = 31;
1489 if (i == 0)
goto S100;
1493 ustar = u - (double)i;
1496 if (ustar <= *(t + i - 1))
goto S60;
1497 w = (ustar - *(t + i - 1)) * *(h + i - 1);
1504 if (s == 1.0) snorm_mt = -y;
1511 w = u * (*(a + i) - aa);
1512 tt = (0.5 * w + aa) * w;
1516 ustar = ranf_mt(tn);
1518 if (ustar > tt)
goto S50;
1520 if (ustar >= u)
goto S70;
1521 ustar = ranf_mt(tn);
1535 if (u < 1.0)
goto S110;
1538 w = u * *(d + i - 1);
1539 tt = (0.5 * w + aa) * w;
1544 ustar = ranf_mt(tn);
1545 if (ustar > tt)
goto S50;
1547 if (ustar >= u)
goto S150;
1551 double fsign(
double num,
double sign)
1554 if ((sign > 0.0f && num < 0.0f) || (sign < 0.0f && num>0.0f))
1563 double gen_norm_mt(
double mu,
double sd,
int tn)
1569 u = 2 * ranf_mt(tn) - 1;
1570 v = 2 * ranf_mt(tn) - 1;
1574 }
while (S >= 1 || S == 0);
1577 x = u * sqrt((-2 * log(S)) / S);
1592 double gen_gamma_mt(
double beta,
double alpha,
int tn)
1594 double d, c, u, v, z, f, alpha2, gamma;
1597 if ((beta <= 0) || (alpha <= 0))
1599 ERR_CRITICAL(
"Gamma distribution parameters in undefined range!\n");
1605 d = alpha - (1.0 / 3.0);
1606 c = 1.0 / (3.0 * sqrt(d));
1611 z = gen_norm_mt(0, 1, tn);
1614 }
while ((z <= (-1.0 / c)) || (log(u) >= (0.5 * z * z + d - d * v + d * log(v))));
1618 return (d * v) / beta;
1629 d = alpha2 - (1.0 / 3.0);
1630 c = 1.0 / (3.0 * sqrt(d));
1635 z = gen_norm_mt(0, 1, tn);
1638 }
while ((z <= (-1.0 / c)) || (log(u) >= (0.5 * z * z + d - d * v + d * log(v))));
1642 gamma = (d * v) / beta;
1649 f = pow(ranf_mt(tn), (1.0 / alpha));
1665 double gen_lognormal(
double mu,
double sigma)
1667 double randnorm, location, scale;
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);
1674 void SampleWithoutReplacement(
int tn,
int k,
int n)
1684 double t, r, a, mu, f;
1689 for (i = 0; i < k; i++)
1693 SamplingQueue[tn][i] = (int)(ranf_mt(tn) * ((double)n));
1701 q = (SamplingQueue[tn][i] == SamplingQueue[tn][j]);
1708 for (i = 0; i < n; i++)
1709 SamplingQueue[tn][i] = i;
1710 for (i = n; i > k; i--)
1712 j = (int)(ranf_mt(tn) * ((double)i));
1715 b = SamplingQueue[tn][j];
1716 SamplingQueue[tn][j] = SamplingQueue[tn][i - 1];
1717 SamplingQueue[tn][i - 1] = b;
1724 for (i = 0; i < n; i++)
1725 SamplingQueue[tn][i] = i;
1726 for (i = 0; i < k; i++)
1728 j = (int)(ranf_mt(tn) * ((double)(n - i)));
1731 b = SamplingQueue[tn][i];
1732 SamplingQueue[tn][i] = SamplingQueue[tn][i + j];
1733 SamplingQueue[tn][i + j] = b;
1743 a = sqrt(log(1 + t / 2 * PI));
1744 a = a + a * a / (3 * r);
1746 b = 2 * MAX_PLACE_SIZE;
1747 f = -1 / (log(1 - mu / ((
double)n)));
1751 i += (int)ceil(-log(ranf_mt(tn)) * f);
1754 SamplingQueue[tn][q] = i - 1;
1756 if (q >= b) i = q = 0;
1803 i = (int)(ranf_mt(tn) * ((double)q));
1804 if (i < q - 1) SamplingQueue[tn][i] = SamplingQueue[tn][q - 1];