+#endif
+}
+
+static void contrast(bcf_p1aux_t *ma, double pc[4]) // mc_cal_y() must be called before hand
+{
+ int k, n1 = ma->n1, n2 = ma->n - ma->n1;
+ long double sum = -1., x, sum_alt;
+ double y;
+ pc[0] = pc[1] = pc[2] = pc[3] = -1.;
+ if (n1 <= 0 || n2 <= 0) return;
+#ifdef _BCF_QUAD
+ { // FIXME: can be improved by skipping zero cells
+ int k1, k2;
+ long double z[3];
+ z[0] = z[1] = z[2] = 0.;
+ for (k1 = 0; k1 <= 2*n1; ++k1)
+ for (k2 = 0; k2 <= 2*n2; ++k2) {
+ double zz = ma->phi[k1+k2] * ma->z1[k1] * ma->z2[k2] * ma->k1k2[k1*(2*n2+1)+k2];
+ if ((double)k1/n1 < (double)k2/n2) z[0] += zz;
+ else if ((double)k1/n1 > (double)k2/n2) z[1] += zz;
+ else z[2] += zz;
+ }
+ sum = z[0] + z[1] + z[2];
+ pc[2] = z[0] / sum; pc[3] = z[1] / sum;
+ }
+#endif
+ for (k = 0, sum_alt = 0.; k <= ma->M; ++k)
+ sum_alt += (long double)ma->phi[k] * ma->z[k];
+// printf("* %lg, %lg *\n", (double)sum, (double)sum_alt); // DEBUG: sum should equal sum_alt
+ sum = sum_alt;
+ y = lgamma(2*n1 + 1) - lgamma(ma->M + 1);
+ for (k = 1, x = 0.; k <= 2 * n2; ++k)
+ x += ma->phi[k] * ma->z2[k] * exp(lgamma(ma->M - k + 1) - lgamma(2*n2 - k + 1) + y);
+ pc[0] = ma->z1[0] * x / sum;
+ y = lgamma(2*n2 + 1) - lgamma(ma->M + 1);
+ for (k = 1, x = 0.; k <= 2 * n1; ++k)
+ x += ma->phi[k] * ma->z1[k] * exp(lgamma(ma->M - k + 1) - lgamma(2*n1 - k + 1) + y);
+ pc[1] = ma->z2[0] * x / sum;
+ for (k = 0; k < 4; ++k) {
+ y = 1. - pc[k];
+ if (y <= 0.) y = 1e-100;
+ pc[k] = (int)(-3.434 * log(y) + .499);
+ if (pc[k] > 99.) pc[k] = 99.;
+ }