+
+double mc_ref_prob(mc_aux_t *ma)
+{
+ int k, i, g;
+ long double PD = 0., Pref = 0.;
+ for (k = 0; k <= ma->n * 2; ++k) {
+ long double x = 1.;
+ double *bk = ma->beta + k * 3;
+ for (i = 0; i < ma->n; ++i) {
+ double y = 0., *pdg = ma->pdg + i * 3;
+ for (g = 0; g < 3; ++g)
+ y += pdg[g] * bk[g];
+ x *= y;
+ }
+ PD += x * ma->alpha[k];
+ }
+ for (k = 0; k <= ma->n * 2; ++k) {
+ long double x = 1.0;
+ for (i = 0; i < ma->n; ++i)
+ x *= ma->pdg[i * 3 + 2] * ma->beta[k * 3 + 2];
+ Pref += x * ma->alpha[k];
+ }
+ return Pref / PD;
+}
+
+int mc_call_gt(const mc_aux_t *ma, double f0, int k)
+{
+ double sum, g[3];
+ double max, f3[3], *pdg = ma->pdg + k * 3;
+ int q, i, max_i;
+ f3[0] = (1.-f0)*(1.-f0); f3[1] = 2.*f0*(1.-f0); f3[2] = f0*f0;
+ for (i = 0, sum = 0.; i < 3; ++i)
+ sum += (g[i] = pdg[i] * f3[i]);
+ for (i = 0, max = -1., max_i = 0; i < 3; ++i) {
+ g[i] /= sum;
+ if (g[i] > max) max = g[i], max_i = i;
+ }
+ max = 1. - max;
+ if (max < 1e-308) max = 1e-308;
+ q = (int)(-3.434 * log(max) + .499);
+ if (q > 99) q = 99;
+ return q<<2|max_i;
+}