pfqfn.cpp

Quantitative Finance C++ source code, asian option pricer

		}
	    }
	}
/* L50: */
    }

/*     SCREENING OF DENOMINATOR ARGUMENTS FOR ZEROES OR NEGATIVE INTEGERS. */

    i__1 = *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {
	if (cr[i1 - 1] == consts_1.zero && cr2[i1 - 1] == consts_1.zero && ci[
		i1 - 1] == consts_1.zero && ci2[i1 - 1] == consts_1.zero) {
	    io___49.ciunit = io_1.nout;
	    s_wsfe(&io___49);
	    do_fio(&c__1, (char *)&i1, (ftnlen)sizeof(integer));
	    e_wsfe();
	    s_stop("", (ftnlen)0);
	}
	i__2 = i1;
	if (ci[i1 - 1] == consts_1.zero && ci2[i1 - 1] == consts_1.zero && b[
		i__2].r < consts_1.zero) {
	    i__2 = i1;
	    i__3 = i1;
	    d__2 = b[i__3].r;
	    i__4 = -nmach;
	    i__5 = i1;
	    d__3 = b[i__5].r;
	    if ((d__1 = b[i__2].r - (doublereal) i_dnnt(&d__2), abs(d__1)) < 
		    pow_di(&consts_1.ten, &i__4) && (icount >= -i_dnnt(&d__3) 
		    || icount == -1)) {
		io___50.ciunit = io_1.nout;
		s_wsfe(&io___50);
		do_fio(&c__1, (char *)&i1, (ftnlen)sizeof(integer));
		e_wsfe();
		s_stop("", (ftnlen)0);
	    }
	}
/* L60: */
    }

    d__1 = pow_di(&consts_1.two, &ibit);
    nmach = (integer) d_lg10(&d__1);
/* Computing MIN */
    d__1 = pow_di(&consts_1.two, &ibit);
    i__1 = *nsigfig, i__2 = (integer) d_lg10(&d__1);
    *nsigfig = min(i__1,i__2);
    i__1 = -(*nsigfig);
    accy = pow_di(&consts_1.ten, &i__1);
    l = ipremax_(&a[1], &b[1], ip, iq, z__);
    if (l == 1) {
	goto L110;
    }

/*     First, estimate the exponent of the maximum term in the pFq series. */

    expon = consts_1.zero;
    xl = (doublereal) l;
    i__1 = *ip;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = i__;
	z__3.r = a[i__2].r + xl, z__3.i = a[i__2].i;
	z__2.r = z__3.r - consts_1.one, z__2.i = z__3.i;
	factor_(&z__1, &z__2);
	i__3 = i__;
	z__5.r = a[i__3].r - consts_1.one, z__5.i = a[i__3].i;
	factor_(&z__4, &z__5);
	expon = expon + z__1.r - z__4.r;
/* L70: */
    }
    i__1 = *iq;
    for (i__ = 1; i__ <= i__1; ++i__) {
	i__2 = i__;
	z__3.r = b[i__2].r + xl, z__3.i = b[i__2].i;
	z__2.r = z__3.r - consts_1.one, z__2.i = z__3.i;
	factor_(&z__1, &z__2);
	i__3 = i__;
	z__5.r = b[i__3].r - consts_1.one, z__5.i = b[i__3].i;
	factor_(&z__4, &z__5);
	expon = expon - z__1.r + z__4.r;
/* L80: */
    }
    z_log(&z__1, z__);
    z__3.r = xl, z__3.i = consts_1.zero;
    factor_(&z__2, &z__3);
    expon = expon + xl * z__1.r - z__2.r;
    d__1 = exp(consts_1.one);
    lmax = (integer) (d_lg10(&d__1) * expon);
    l = lmax;

/*     Now, estimate the exponent of where the pFq series will terminate. */

    z__1.r = consts_1.one, z__1.i = consts_1.zero;
    temp1.r = z__1.r, temp1.i = z__1.i;
    creal = consts_1.one;
    i__1 = *ip;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = i1 - 1;
	i__3 = i1 - 1;
	z__3.r = ar[i__2], z__3.i = ai[i__3];
	z__2.r = temp1.r * z__3.r - temp1.i * z__3.i, z__2.i = temp1.r * 
		z__3.i + temp1.i * z__3.r;
	z__1.r = z__2.r / sigfig, z__1.i = z__2.i / sigfig;
	temp1.r = z__1.r, temp1.i = z__1.i;
/* L90: */
    }
    i__1 = *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = i1 - 1;
	i__3 = i1 - 1;
	z__3.r = cr[i__2], z__3.i = ci[i__3];
	z__2.r = z__3.r / sigfig, z__2.i = z__3.i / sigfig;
	z_div(&z__1, &temp1, &z__2);
	temp1.r = z__1.r, temp1.i = z__1.i;
	creal *= cr[i1 - 1];
/* L100: */
    }
    z__2.r = xr, z__2.i = xi;
    z__1.r = temp1.r * z__2.r - temp1.i * z__2.i, z__1.i = temp1.r * z__2.i + 
	    temp1.i * z__2.r;
    temp1.r = z__1.r, temp1.i = z__1.i;

/*     Triple it to make sure. */

    l *= 3;

/*     Divide the number of significant figures necessary by the number of */
/*     digits available per array position. */

    l = ((l << 1) + *nsigfig) / nmach + 2;

/*     Make sure there are at least 5 array positions used. */

L110:
    l = max(l,5);
    l = max(l,*ix);
/*     write (6,*) ' Estimated value of L=',L */
    if (l < 0 || l > 25000) {
	io___59.ciunit = io_1.nout;
	s_wsfe(&io___59);
	do_fio(&c__1, (char *)&c__25000, (ftnlen)sizeof(integer));
	e_wsfe();
	s_stop("", (ftnlen)0);
    }
    if (*nsigfig > nmach) {
	io___60.ciunit = io_1.nout;
	s_wsfe(&io___60);
	do_fio(&c__1, (char *)&nmach, (ftnlen)sizeof(integer));
	e_wsfe();
    }
    sumr[0] = consts_1.one;
    sumi[0] = consts_1.one;
    numr[0] = consts_1.one;
    numi[0] = consts_1.one;
    denomr[0] = consts_1.one;
    denomi[0] = consts_1.one;
    i__1 = l + 1;
    for (i__ = 0; i__ <= i__1; ++i__) {
	sumr[i__ + 1] = consts_1.zero;
	sumi[i__ + 1] = consts_1.zero;
	numr[i__ + 1] = consts_1.zero;
	numi[i__ + 1] = consts_1.zero;
	denomr[i__ + 1] = consts_1.zero;
	denomi[i__ + 1] = consts_1.zero;
/* L120: */
    }
    sumr[2] = consts_1.one;
    numr[2] = consts_1.one;
    denomr[2] = consts_1.one;
    cnt = sigfig;
    z__1.r = consts_1.zero, z__1.i = consts_1.zero;
    temp.r = z__1.r, temp.i = z__1.i;
    oldtemp.r = temp.r, oldtemp.i = temp.i;
    ixcnt = 0;
    rexp = ibit / 2;
    x = rexp * (sumr[l + 2] - 2);
    rr10 = x * log2;
    ir10 = (integer) rr10;
    rr10 -= ir10;
    x = rexp * (sumi[l + 2] - 2);
    ri10 = x * log2;
    ii10 = (integer) ri10;
    ri10 -= ii10;
    d__1 = sumr[2] * rmax * rmax + sumr[3] * rmax + sumr[4];
    dum1 = d_sign(&d__1, sumr);
    d__1 = sumi[2] * rmax * rmax + sumi[3] * rmax + sumi[4];
    dum2 = d_sign(&d__1, sumi);
    dum1 *= pow_dd(&c_b64, &rr10);
    dum2 *= pow_dd(&c_b64, &ri10);
    z__1.r = dum1, z__1.i = dum2;
    cdum1.r = z__1.r, cdum1.i = z__1.i;
    x = rexp * (denomr[l + 2] - 2);
    rr10 = x * log2;
    ir10 = (integer) rr10;
    rr10 -= ir10;
    x = rexp * (denomi[l + 2] - 2);
    ri10 = x * log2;
    ii10 = (integer) ri10;
    ri10 -= ii10;
    d__1 = denomr[2] * rmax * rmax + denomr[3] * rmax + denomr[4];
    dum1 = d_sign(&d__1, denomr);
    d__1 = denomi[2] * rmax * rmax + denomi[3] * rmax + denomi[4];
    dum2 = d_sign(&d__1, denomi);
    dum1 *= pow_dd(&c_b64, &rr10);
    dum2 *= pow_dd(&c_b64, &ri10);
    z__1.r = dum1, z__1.i = dum2;
    cdum2.r = z__1.r, cdum2.i = z__1.i;
    z_div(&z__1, &cdum1, &cdum2);
    temp.r = z__1.r, temp.i = z__1.i;
/*     130 IF (IP .GT. 0) THEN */
L130:
    if (*ip < 0) {
	if (sumr[2] < consts_1.half) {
	    mx1 = sumi[l + 2];
	} else if (sumi[2] < consts_1.half) {
	    mx1 = sumr[l + 2];
	} else {
/* Computing MAX */
	    d__1 = sumr[l + 2], d__2 = sumi[l + 2];
	    mx1 = max(d__1,d__2);
	}
	if (numr[2] < consts_1.half) {
	    mx2 = numi[l + 2];
	} else if (numi[2] < consts_1.half) {
	    mx2 = numr[l + 2];
	} else {
/* Computing MAX */
	    d__1 = numr[l + 2], d__2 = numi[l + 2];
	    mx2 = max(d__1,d__2);
	}
	if (mx1 - mx2 > 2.f) {
	    if (creal >= consts_1.zero) {
/*     write (6,*) ' cdabs(temp1/cnt)=',cdabs(temp1/cnt) */
		z__1.r = temp1.r / cnt, z__1.i = temp1.i / cnt;
		if (z_abs(&z__1) <= consts_1.one) {
		    goto L240;
		}
	    }
	}
    } else {
	arydiv_(sumr, sumi, denomr, denomi, &temp, &l, lnpfq, &rmax, &ibit);

/*     First, estimate the exponent of the maximum term in the pFq series. */

	expon = consts_1.zero;
	xl = (doublereal) ixcnt;
	i__1 = *ip;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    i__2 = i__;
	    z__3.r = a[i__2].r + xl, z__3.i = a[i__2].i;
	    z__2.r = z__3.r - consts_1.one, z__2.i = z__3.i;
	    factor_(&z__1, &z__2);
	    i__3 = i__;
	    z__5.r = a[i__3].r - consts_1.one, z__5.i = a[i__3].i;
	    factor_(&z__4, &z__5);
	    expon = expon + z__1.r - z__4.r;
/* L140: */
	}
	i__1 = *iq;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    i__2 = i__;
	    z__3.r = b[i__2].r + xl, z__3.i = b[i__2].i;
	    z__2.r = z__3.r - consts_1.one, z__2.i = z__3.i;
	    factor_(&z__1, &z__2);
	    i__3 = i__;
	    z__5.r = b[i__3].r - consts_1.one, z__5.i = b[i__3].i;
	    factor_(&z__4, &z__5);
	    expon = expon - z__1.r + z__4.r;
/* L150: */
	}
	z_log(&z__1, z__);
	z__3.r = xl, z__3.i = consts_1.zero;
	factor_(&z__2, &z__3);
	expon = expon + xl * z__1.r - z__2.r;
	d__1 = exp(consts_1.one);
	lmax = (integer) (d_lg10(&d__1) * expon);
	z__1.r = oldtemp.r - temp.r, z__1.i = oldtemp.i - temp.i;
	z__2.r = accy * temp.r, z__2.i = accy * temp.i;
	if (z_abs(&z__1) < z_abs(&z__2)) {
	    goto L240;
	}
	oldtemp.r = temp.r, oldtemp.i = temp.i;
    }
    if (ixcnt == icount) {
	goto L240;
    }
    ++ixcnt;
    i__1 = *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {

/*     TAKE THE CURRENT SUM AND MULTIPLY BY THE DENOMINATOR OF THE NEXT */
/*     TERM, FOR BOTH THE MOST SIGNIFICANT HALF (CR,CI) AND THE LEAST */
/*     SIGNIFICANT HALF (CR2,CI2). */

	cmpmul_(sumr, sumi, &cr[i1 - 1], &ci[i1 - 1], qr1, qi1, wk1, wk2, wk3,
		 wk4, wk5, wk6, &l, &rmax);
	cmpmul_(sumr, sumi, &cr2[i1 - 1], &ci2[i1 - 1], qr2, qi2, wk1, wk2, 
		wk3, wk4, wk5, wk6, &l, &rmax);
	qr2[l + 2] += -1;
	qi2[l + 2] += -1;

/*     STORE THIS TEMPORARILY IN THE SUM ARRAYS. */

	cmpadd_(qr1, qi1, qr2, qi2, sumr, sumi, wk1, &l, &rmax);
/* L160: */
    }

/*     MULTIPLY BY THE FACTORIAL TERM. */

    armult_(sumr, &cnt, sumr, wk6, &l, &rmax);
    armult_(sumi, &cnt, sumi, wk6, &l, &rmax);

/*     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT. */

    i__1 = *ip - *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {
	armult_(sumr, &sigfig, sumr, wk6, &l, &rmax);
	armult_(sumi, &sigfig, sumi, wk6, &l, &rmax);
/* L170: */
    }
    i__1 = *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {

/*     UPDATE THE DENOMINATOR. */

	cmpmul_(denomr, denomi, &cr[i1 - 1], &ci[i1 - 1], qr1, qi1, wk1, wk2, 
		wk3, wk4, wk5, wk6, &l, &rmax);
	cmpmul_(denomr, denomi, &cr2[i1 - 1], &ci2[i1 - 1], qr2, qi2, wk1, 
		wk2, wk3, wk4, wk5, wk6, &l, &rmax);
	qr2[l + 2] += -1;
	qi2[l + 2] += -1;
	cmpadd_(qr1, qi1, qr2, qi2, denomr, denomi, wk1, &l, &rmax);
/* L180: */
    }

/*     MULTIPLY BY THE FACTORIAL TERM. */

    armult_(denomr, &cnt, denomr, wk6, &l, &rmax);
    armult_(denomi, &cnt, denomi, wk6, &l, &rmax);

/*     MULTIPLY BY THE SCALING FACTOR, SIGFIG, TO KEEP THE SCALE CORRECT. */

    i__1 = *ip - *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {
	armult_(denomr, &sigfig, denomr, wk6, &l, &rmax);
	armult_(denomi, &sigfig, denomi, wk6, &l, &rmax);
/* L190: */
    }

/*     FORM THE NEXT NUMERATOR TERM BY MULTIPLYING THE CURRENT */
/*     NUMERATOR TERM (AN ARRAY) WITH THE A ARGUMENT (A SCALAR). */

    i__1 = *ip;
    for (i1 = 1; i1 <= i__1; ++i1) {
	cmpmul_(numr, numi, &ar[i1 - 1], &ai[i1 - 1], qr1, qi1, wk1, wk2, wk3,
		 wk4, wk5, wk6, &l, &rmax);
	cmpmul_(numr, numi, &ar2[i1 - 1], &ai2[i1 - 1], qr2, qi2, wk1, wk2, 
		wk3, wk4, wk5, wk6, &l, &rmax);
	qr2[l + 2] += -1;
	qi2[l + 2] += -1;
	cmpadd_(qr1, qi1, qr2, qi2, numr, numi, wk1, &l, &rmax);
/* L200: */
    }

/*     FINISH THE NEW NUMERATOR TERM BY MULTIPLYING BY THE Z ARGUMENT. */

    cmpmul_(numr, numi, &xr, &xi, qr1, qi1, wk1, wk2, wk3, wk4, wk5, wk6, &l,