-/*
- if (ma->n1 > 0 && ma->n1 < ma->n) {
- int i;
- double y[5];
- printf("*****\n");
- for (i = 0; i <= 2; ++i)
- printf("(%lf,%lf) ", ma->z1[i], ma->z2[i]);
- printf("\n");
- y[0] = ma->z1[0] * ma->z2[0];
- y[1] = 1./2. * (ma->z1[0] * ma->z2[1] + ma->z1[1] * ma->z2[0]);
- y[2] = 1./6. * (ma->z1[0] * ma->z2[2] + ma->z1[2] * ma->z2[0]) + 4./6. * ma->z1[1] * ma->z2[1];
- y[3] = 1./2. * (ma->z1[1] * ma->z2[2] + ma->z1[2] * ma->z2[1]);
- y[4] = ma->z1[2] * ma->z2[2];
- for (i = 0; i <= 4; ++i) printf("(%lf,%lf) ", ma->z[i], y[i]);
- printf("\n");
+}
+
+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 sum1, sum2;
+ pc[0] = pc[1] = pc[2] = pc[3] = -1.;
+ if (n1 <= 0 || n2 <= 0) return;
+ 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;