+static void mc_cal_y(bcf_p1aux_t *ma)
+{
+ if (ma->n1 > 0 && ma->n1 < ma->n) {
+ int k;
+ double x;
+ memset(ma->z1, 0, sizeof(double) * (2 * ma->n1 + 1));
+ memset(ma->z2, 0, sizeof(double) * (2 * (ma->n - ma->n1) + 1));
+ ma->t1 = ma->t2 = 0.;
+ mc_cal_y_core(ma, ma->n1);
+ ma->t2 = ma->t;
+ memcpy(ma->z2, ma->z, sizeof(double) * (2 * (ma->n - ma->n1) + 1));
+ mc_cal_y_core(ma, 0);
+ // rescale z
+ x = exp(ma->t - (ma->t1 + ma->t2));
+ for (k = 0; k <= ma->M; ++k) ma->z[k] *= x;
+ } else mc_cal_y_core(ma, 0);
+/*
+ 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");
+ }
+*/
+}
+