-#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.;
+ for (k = 0, sum1 = 0.; k <= 2*n1; ++k) sum1 += ma->phi1[k] * ma->z1[k];
+ for (k = 0, sum2 = 0.; k <= 2*n2; ++k) sum2 += ma->phi2[k] * ma->z2[k];
+ pc[2] = ma->phi1[2*n1] * ma->z1[2*n1] / sum1;
+ pc[3] = ma->phi2[2*n2] * ma->z2[2*n2] / sum2;
+ for (k = 2; k < 4; ++k) {
+ pc[k] = pc[k] > .5? -(-4.343 * log(1. - pc[k] + TINY) + .499) : -4.343 * log(pc[k] + TINY) + .499;
+ pc[k] = (int)pc[k];
+ if (pc[k] > 99) pc[k] = 99;
+ if (pc[k] < -99) pc[k] = -99;