-#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;
- }
-#else
- pc[2] = pc[3] = 0.;
-#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;
- // the variant is specific to group2
-// printf("%lg %lg %lg %lg\n", ma->z[2*(n1+n2)]/exp(ma->t - (ma->t1 + ma->t2)), ma->z1[2*n1], ma->z2[2*n2], (double)sum);
- y = lgamma(2*n2 + 1) - lgamma(ma->M + 1);
- for (k = 0, x = 0.; k < 2 * n2; ++k)
- x += ma->phi[2*n1+k] * ma->z2[k] * expl(lgamma(2*n1 + k + 1) - lgamma(k + 1) + y);
- pc[1] = ma->z1[2*n1] * x / sum;
- for (k = 1, x = 0.; k <= 2 * n2; ++k)
- x += ma->phi[k] * ma->z2[k] * expl(lgamma(ma->M - k + 1) - lgamma(2*n2 - k + 1) + y);
- pc[1] += ma->z1[0] * x / sum;
- // the variant is specific to group1
- y = lgamma(2*n1 + 1) - lgamma(ma->M + 1);
- for (k = 0, x = 0.; k < 2 * n1; ++k)
- x += ma->phi[2*n2+k] * ma->z1[k] * expl(lgamma(2*n2 + k + 1) - lgamma(k + 1) + y);
- pc[0] = ma->z2[2*n2] * x / sum;
- for (k = 1, x = 0.; k <= 2 * n1; ++k)
- x += ma->phi[k] * ma->z1[k] * expl(lgamma(ma->M - k + 1) - lgamma(2*n1 - k + 1) + y);
- pc[0] += ma->z2[0] * x / sum;
- // rescale