pfqfn.cpp

Quantitative Finance C++ source code, asian option pricer

	    &rmax);
    cmpmul_(numr, numi, &xr2, &xi2, 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);

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

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

/*     FINALLY, ADD THE NEW NUMERATOR TERM WITH THE CURRENT RUNNING */
/*     SUM OF THE NUMERATOR AND STORE THE NEW RUNNING SUM IN SUMR, SUMI. */

    cmpadd_(sumr, sumi, numr, numi, sumr, sumi, wk1, &l, &rmax);

/*     BECAUSE SIGFIG REPRESENTS "ONE" ON THE NEW SCALE, ADD SIGFIG */
/*     TO THE CURRENT COUNT AND, CONSEQUENTLY, TO THE IP ARGUMENTS */
/*     IN THE NUMERATOR AND THE IQ ARGUMENTS IN THE DENOMINATOR. */

    cnt += sigfig;
    i__1 = *ip;
    for (i1 = 1; i1 <= i__1; ++i1) {
	ar[i1 - 1] += sigfig;
/* L220: */
    }
    i__1 = *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {
	cr[i1 - 1] += sigfig;
/* L230: */
    }
    goto L130;
L240:
    arydiv_(sumr, sumi, denomr, denomi, &final, &l, lnpfq, &rmax, &ibit);
/*     write (6,*) 'Number of terms=',ixcnt */
     ret_val->r = final.r,  ret_val->i = final.i;
    return ;
} /* hyper_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ARADD                             * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Accepts two arrays of numbers and returns     * */
/*     *    the sum of the array.  Each array is holding the value    * */
/*     *    of one number in the series.  The parameter L is the      * */
/*     *    size of the array representing the number and RMAX is     * */
/*     *    the actual number of digits needed to give the numbers    * */
/*     *    the desired accuracy.                                     * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int aradd_(doublereal *a, doublereal *b, doublereal *c__, 
	doublereal *z__, integer *l, doublereal *rmax)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1;

    /* Builtin functions */
    double d_nint(doublereal *);

    /* Local variables */
    static integer i__, j, ediff;

    /* Parameter adjustments */
    ++z__;
    ++c__;
    ++b;
    ++a;

    /* Function Body */
    i__1 = *l + 1;
    for (i__ = 0; i__ <= i__1; ++i__) {
	z__[i__] = consts_1.zero;
/* L10: */
    }
    d__1 = a[*l + 1] - b[*l + 1];
    ediff = (integer) d_nint(&d__1);
    if (abs(a[1]) < consts_1.half || ediff <= -(*l)) {
	goto L20;
    }
    if (abs(b[1]) < consts_1.half || ediff >= *l) {
	goto L40;
    }
    goto L60;
L20:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
	c__[i__] = b[i__];
/* L30: */
    }
    goto L350;
L40:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
	c__[i__] = a[i__];
/* L50: */
    }
    goto L350;
L60:
    z__[-1] = a[-1];
    if ((d__1 = a[-1] - b[-1], abs(d__1)) < consts_1.half) {
	goto L80;
    }
    if (ediff > 0) {
	z__[*l + 1] = a[*l + 1];
	goto L190;
    }
    if (ediff < 0) {
	z__[*l + 1] = b[*l + 1];
	z__[-1] = b[-1];
	goto L240;
    }
    i__1 = *l;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if (a[i__] > b[i__]) {
	    z__[*l + 1] = a[*l + 1];
	    goto L190;
	}
	if (a[i__] < b[i__]) {
	    z__[*l + 1] = b[*l + 1];
	    z__[-1] = b[-1];
	    goto L240;
	}
/* L70: */
    }
    goto L300;
L80:
    if (ediff > 0) {
	goto L110;
    }
    if (ediff < 0) {
	goto L150;
    }
    z__[*l + 1] = a[*l + 1];
    for (i__ = *l; i__ >= 1; --i__) {
	z__[i__] = a[i__] + b[i__] + z__[i__];
	if (z__[i__] >= *rmax) {
	    z__[i__] -= *rmax;
	    z__[i__ - 1] = consts_1.one;
	}
/* L90: */
    }
    if (z__[0] > consts_1.half) {
	for (i__ = *l; i__ >= 1; --i__) {
	    z__[i__] = z__[i__ - 1];
/* L100: */
	}
	z__[*l + 1] += consts_1.one;
	z__[0] = consts_1.zero;
    }
    goto L300;
L110:
    z__[*l + 1] = a[*l + 1];
    i__1 = ediff + 1;
    for (i__ = *l; i__ >= i__1; --i__) {
	z__[i__] = a[i__] + b[i__ - ediff] + z__[i__];
	if (z__[i__] >= *rmax) {
	    z__[i__] -= *rmax;
	    z__[i__ - 1] = consts_1.one;
	}
/* L120: */
    }
    for (i__ = ediff; i__ >= 1; --i__) {
	z__[i__] = a[i__] + z__[i__];
	if (z__[i__] >= *rmax) {
	    z__[i__] -= *rmax;
	    z__[i__ - 1] = consts_1.one;
	}
/* L130: */
    }
    if (z__[0] > consts_1.half) {
	for (i__ = *l; i__ >= 1; --i__) {
	    z__[i__] = z__[i__ - 1];
/* L140: */
	}
	++z__[*l + 1];
	z__[0] = consts_1.zero;
    }
    goto L300;
L150:
    z__[*l + 1] = b[*l + 1];
    i__1 = 1 - ediff;
    for (i__ = *l; i__ >= i__1; --i__) {
	z__[i__] = a[i__ + ediff] + b[i__] + z__[i__];
	if (z__[i__] >= *rmax) {
	    z__[i__] -= *rmax;
	    z__[i__ - 1] = consts_1.one;
	}
/* L160: */
    }
    for (i__ = -ediff; i__ >= 1; --i__) {
	z__[i__] = b[i__] + z__[i__];
	if (z__[i__] >= *rmax) {
	    z__[i__] -= *rmax;
	    z__[i__ - 1] = consts_1.one;
	}
/* L170: */
    }
    if (z__[0] > consts_1.half) {
	for (i__ = *l; i__ >= 1; --i__) {
	    z__[i__] = z__[i__ - 1];
/* L180: */
	}
	z__[*l + 1] += consts_1.one;
	z__[0] = consts_1.zero;
    }
    goto L300;
L190:
    if (ediff > 0) {
	goto L210;
    }
    for (i__ = *l; i__ >= 1; --i__) {
	z__[i__] = a[i__] - b[i__] + z__[i__];
	if (z__[i__] < consts_1.zero) {
	    z__[i__] += *rmax;
	    z__[i__ - 1] = -consts_1.one;
	}
/* L200: */
    }
    goto L290;
L210:
    i__1 = ediff + 1;
    for (i__ = *l; i__ >= i__1; --i__) {
	z__[i__] = a[i__] - b[i__ - ediff] + z__[i__];
	if (z__[i__] < consts_1.zero) {
	    z__[i__] += *rmax;
	    z__[i__ - 1] = -consts_1.one;
	}
/* L220: */
    }
    for (i__ = ediff; i__ >= 1; --i__) {
	z__[i__] = a[i__] + z__[i__];
	if (z__[i__] < consts_1.zero) {
	    z__[i__] += *rmax;
	    z__[i__ - 1] = -consts_1.one;
	}
/* L230: */
    }
    goto L290;
L240:
    if (ediff < 0) {
	goto L260;
    }
    for (i__ = *l; i__ >= 1; --i__) {
	z__[i__] = b[i__] - a[i__] + z__[i__];
	if (z__[i__] < consts_1.zero) {
	    z__[i__] += *rmax;
	    z__[i__ - 1] = -consts_1.one;
	}
/* L250: */
    }
    goto L290;
L260:
    i__1 = 1 - ediff;
    for (i__ = *l; i__ >= i__1; --i__) {
	z__[i__] = b[i__] - a[i__ + ediff] + z__[i__];
	if (z__[i__] < consts_1.zero) {
	    z__[i__] += *rmax;
	    z__[i__ - 1] = -consts_1.one;
	}
/* L270: */
    }
    for (i__ = -ediff; i__ >= 1; --i__) {
	z__[i__] = b[i__] + z__[i__];
	if (z__[i__] < consts_1.zero) {
	    z__[i__] += *rmax;
	    z__[i__ - 1] = -consts_1.one;
	}
/* L280: */
    }
L290:
    if (z__[1] > consts_1.half) {
	goto L330;
    }
    i__ = 1;
L300:
    ++i__;
    if (z__[i__] < consts_1.half && i__ < *l + 1) {
	goto L300;
    }
    if (i__ == *l + 1) {
	z__[-1] = consts_1.one;
	z__[*l + 1] = consts_1.zero;
	goto L330;
    }
    i__1 = *l + 1 - i__;
    for (j = 1; j <= i__1; ++j) {
	z__[j] = z__[j + i__ - 1];
/* L310: */
    }
    i__1 = *l;
    for (j = *l + 2 - i__; j <= i__1; ++j) {
	z__[j] = consts_1.zero;
/* L320: */
    }
    z__[*l + 1] = z__[*l + 1] - i__ + 1;
L330:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
	c__[i__] = z__[i__];
/* L340: */
    }
L350:
    if (c__[1] < consts_1.half) {
	c__[-1] = consts_1.one;
	c__[*l + 1] = consts_1.zero;
    }
    return 0;
} /* aradd_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ARSUB                             * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Accepts two arrays and subtracts each element * */
/*     *    in the second array from the element in the first array   * */
/*     *    and returns the solution.  The parameters L and RMAX are  * */
/*     *    the size of the array and the number of digits needed for * */
/*     *    the accuracy, respectively.                               * */
/*     *                                                              * */
/*     *  Subprograms called: ARADD                                   * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int arsub_(doublereal *a, doublereal *b, doublereal *c__, 
	doublereal *wk1, doublereal *wk2, integer *l, doublereal *rmax)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    static integer i__;
    extern /* Subroutine */ int aradd_(doublereal *, doublereal *, doublereal 
	    *, doublereal *, integer *, doublereal *);

    /* Parameter adjustments */
    ++wk2;
    ++wk1;
    ++c__;
    ++b;
    ++a;

    /* Function Body */
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
	wk2[i__] = b[i__];
/* L10: */
    }
    wk2[-1] = -consts_1.one * wk2[-1];
    aradd_(&a[-1], &wk2[-1], &c__[-1], &wk1[-1], l, rmax);
    return 0;
} /* arsub_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ARMULT                            * */
/*     *                                                              * */