pfqfn.cpp

Quantitative Finance C++ source code, asian option pricer

/*     *                                                              * */
/*     *  Description : Accepts two arrays and returns the product.   * */
/*     *    L and RMAX are the size of the arrays and the number of   * */
/*     *    digits needed to represent the numbers with the required  * */
/*     *    accuracy.                                                 * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int armult_(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_sign(doublereal *, doublereal *), d_int(doublereal *);

    /* Local variables */
    static integer i__;
    static doublereal b2, carry;

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

    /* Function Body */
    z__[-1] = d_sign(&consts_1.one, b) * a[-1];
    b2 = abs(*b);
    z__[*l + 1] = a[*l + 1];
    i__1 = *l;
    for (i__ = 0; i__ <= i__1; ++i__) {
	z__[i__] = consts_1.zero;
/* L10: */
    }
    if (b2 <= consts_1.eps || a[1] <= consts_1.eps) {
	z__[-1] = consts_1.one;
	z__[*l + 1] = consts_1.zero;
	goto L60;
    }
    for (i__ = *l; i__ >= 1; --i__) {
	z__[i__] = a[i__] * b2 + z__[i__];
	if (z__[i__] >= *rmax) {
	    d__1 = z__[i__] / *rmax;
	    carry = d_int(&d__1);
	    z__[i__] -= carry * *rmax;
	    z__[i__ - 1] = carry;
	}
/* L20: */
    }
    if (z__[0] < consts_1.half) {
	goto L50;
    }
    for (i__ = *l; i__ >= 1; --i__) {
	z__[i__] = z__[i__ - 1];
/* L30: */
    }
    z__[*l + 1] += consts_1.one;
    if (z__[1] >= *rmax) {
	for (i__ = *l; i__ >= 1; --i__) {
	    z__[i__] = z__[i__ - 1];
/* L40: */
	}
	d__1 = z__[1] / *rmax;
	carry = d_int(&d__1);
	z__[2] -= carry * *rmax;
	z__[1] = carry;
	z__[*l + 1] += consts_1.one;
    }
    z__[0] = consts_1.zero;
L50:
L60:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
	c__[i__] = z__[i__];
/* L70: */
    }
    if (c__[1] < consts_1.half) {
	c__[-1] = consts_1.one;
	c__[*l + 1] = consts_1.zero;
    }
    return 0;
} /* armult_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE CMPADD                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Takes two arrays representing one real and    * */
/*     *    one imaginary part, and adds two arrays representing      * */
/*     *    another complex number and returns two array holding the  * */
/*     *    complex sum.                                              * */
/*     *              (CR,CI) = (AR+BR, AI+BI)                        * */
/*     *                                                              * */
/*     *  Subprograms called: ARADD                                   * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int cmpadd_(doublereal *ar, doublereal *ai, doublereal *br, 
	doublereal *bi, doublereal *cr, doublereal *ci, doublereal *wk1, 
	integer *l, doublereal *rmax)
{
    extern /* Subroutine */ int aradd_(doublereal *, doublereal *, doublereal 
	    *, doublereal *, integer *, doublereal *);

    /* Parameter adjustments */
    ++wk1;
    ++ci;
    ++cr;
    ++bi;
    ++br;
    ++ai;
    ++ar;

    /* Function Body */
    aradd_(&ar[-1], &br[-1], &cr[-1], &wk1[-1], l, rmax);
    aradd_(&ai[-1], &bi[-1], &ci[-1], &wk1[-1], l, rmax);
    return 0;
} /* cmpadd_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE CMPSUB                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Takes two arrays representing one real and    * */
/*     *    one imaginary part, and subtracts two arrays representing * */
/*     *    another complex number and returns two array holding the  * */
/*     *    complex sum.                                              * */
/*     *              (CR,CI) = (AR+BR, AI+BI)                        * */
/*     *                                                              * */
/*     *  Subprograms called: ARADD                                   * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int cmpsub_(doublereal *ar, doublereal *ai, doublereal *br, 
	doublereal *bi, doublereal *cr, doublereal *ci, doublereal *wk1, 
	doublereal *wk2, integer *l, doublereal *rmax)
{
    extern /* Subroutine */ int arsub_(doublereal *, doublereal *, doublereal 
	    *, doublereal *, doublereal *, integer *, doublereal *);

    /* Parameter adjustments */
    ++wk2;
    ++wk1;
    ++ci;
    ++cr;
    ++bi;
    ++br;
    ++ai;
    ++ar;

    /* Function Body */
    arsub_(&ar[-1], &br[-1], &cr[-1], &wk1[-1], &wk2[-1], l, rmax);
    arsub_(&ai[-1], &bi[-1], &ci[-1], &wk1[-1], &wk2[-1], l, rmax);
    return 0;
} /* cmpsub_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE CMPMUL                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Takes two arrays representing one real and    * */
/*     *    one imaginary part, and multiplies it with two arrays     * */
/*     *    representing another complex number and returns the       * */
/*     *    complex product.                                          * */
/*     *                                                              * */
/*     *  Subprograms called: ARMULT, ARSUB, ARADD                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int cmpmul_(doublereal *ar, doublereal *ai, doublereal *br, 
	doublereal *bi, doublereal *cr, doublereal *ci, doublereal *wk1, 
	doublereal *wk2, doublereal *cr2, doublereal *d1, doublereal *d2, 
	doublereal *wk6, 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 *), arsub_(doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, integer *,
	     doublereal *), armult_(doublereal *, doublereal *, doublereal *, 
	    doublereal *, integer *, doublereal *);

    /* Parameter adjustments */
    ++d2;
    ++d1;
    ++cr2;
    ++wk2;
    ++wk1;
    ++ci;
    ++cr;
    ++ai;
    ++ar;

    /* Function Body */
    armult_(&ar[-1], br, &d1[-1], wk6, l, rmax);
    armult_(&ai[-1], bi, &d2[-1], wk6, l, rmax);
    arsub_(&d1[-1], &d2[-1], &cr2[-1], &wk1[-1], &wk2[-1], l, rmax);
    armult_(&ar[-1], bi, &d1[-1], wk6, l, rmax);
    armult_(&ai[-1], br, &d2[-1], wk6, l, rmax);
    aradd_(&d1[-1], &d2[-1], &ci[-1], &wk1[-1], l, rmax);
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
	cr[i__] = cr2[i__];
/* L10: */
    }
    return 0;
} /* cmpmul_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ARYDIV                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Returns the double precision complex number   * */
/*     *    resulting from the division of four arrays, representing  * */
/*     *    two complex numbers.  The number returned will be in one  * */
/*     *    of two different forms:  either standard scientific or as * */
/*     *    the log (base 10) of the number.                          * */
/*     *                                                              * */
/*     *  Subprograms called: CONV21, CONV12, EADD, ECPDIV, EMULT.    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int arydiv_(doublereal *ar, doublereal *ai, doublereal *br, 
	doublereal *bi, doublecomplex *c__, integer *l, integer *lnpfq, 
	doublereal *rmax, integer *ibit)
{
    /* System generated locals */
    doublereal d__1;
    doublecomplex z__1;

    /* Builtin functions */
    //double d_lg10(doublereal *), d_sign(doublereal *, doublereal *), pow_dd(
	  //  doublereal *, doublereal *), atan2(doublereal, doublereal), log(
	  //  doublereal);

    /* Local variables */
    static doublereal x, e1, e2, e3, n1, n2, n3, x1, x2, ae[4]	/* was [2][2] 
	    */, be[4]	/* was [2][2] */, ce[4]	/* was [2][2] */;
    static integer ii10, ir10;
    static doublereal ri10, phi, rr10, dum1, dum2;
    extern /* Subroutine */ int eadd_(doublereal *, doublereal *, doublereal *
	    , doublereal *, doublereal *, doublereal *);
    static doublecomplex cdum;
    static real dnum;
    static integer rexp;
    extern /* Subroutine */ int conv12_(doublecomplex *, doublereal *), 
	    conv21_(doublereal *, doublecomplex *), emult_(doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *), ecpdiv_(doublereal *, doublereal *, doublereal *);
    static doublereal tenmax;
    static integer itnmax;

    /* Parameter adjustments */
    ++bi;
    ++br;
    ++ai;
    ++ar;

    /* Function Body */
    rexp = *ibit / 2;
    x = rexp * (ar[*l + 1] - 2);
    rr10 = x * d_lg10(&consts_1.two) / d_lg10(&consts_1.ten);
    ir10 = (integer) rr10;
    rr10 -= ir10;
    x = rexp * (ai[*l + 1] - 2);
    ri10 = x * d_lg10(&consts_1.two) / d_lg10(&consts_1.ten);
    ii10 = (integer) ri10;
    ri10 -= ii10;
    d__1 = ar[1] * *rmax * *rmax + ar[2] * *rmax + ar[3];
    dum1 = d_sign(&d__1, &ar[-1]);
    d__1 = ai[1] * *rmax * *rmax + ai[2] * *rmax + ai[3];
    dum2 = d_sign(&d__1, &ai[-1]);
    dum1 *= pow_dd(&c_b64, &rr10);
    dum2 *= pow_dd(&c_b64, &ri10);
    z__1.r = dum1, z__1.i = dum2;
    cdum.r = z__1.r, cdum.i = z__1.i;
    conv12_(&cdum, ae);
    ae[2] += ir10;
    ae[3] += ii10;
    x = rexp * (br[*l + 1] - 2);
    rr10 = x * d_lg10(&consts_1.two) / d_lg10(&consts_1.ten);
    ir10 = (integer) rr10;
    rr10 -= ir10;
    x = rexp * (bi[*l + 1] - 2);
    ri10 = x * d_lg10(&consts_1.two) / d_lg10(&consts_1.ten);
    ii10 = (integer) ri10;
    ri10 -= ii10;
    d__1 = br[1] * *rmax * *rmax + br[2] * *rmax + br[3];
    dum1 = d_sign(&d__1, &br[-1]);
    d__1 = bi[1] * *rmax * *rmax + bi[2] * *rmax + bi[3];
    dum2 = d_sign(&d__1, &bi[-1]);
    dum1 *= pow_dd(&c_b64, &rr10);
    dum2 *= pow_dd(&c_b64, &ri10);
    z__1.r = dum1, z__1.i = dum2;
    cdum.r = z__1.r, cdum.i = z__1.i;
    conv12_(&cdum, be);
    be[2] += ir10;
    be[3] += ii10;
    ecpdiv_(ae, be, ce);
    if (*lnpfq == 0) {
	conv21_(ce, c__);
    } else {
	emult_(ce, &ce[2], ce, &ce[2], &n1, &e1);
	emult_(&ce[1], &ce[3], &ce[1], &ce[3], &n2, &e2);
	eadd_(&n1, &e1, &n2, &e2, &n3, &e3);
	n1 = ce[0];
	e1 = ce[2] - ce[3];
	x2 = ce[1];

/*     TENMAX - MAXIMUM SIZE OF EXPONENT OF 10 */
/*     THE FOLLOWING CODE CAN BE USED TO DETERMINE TENMAX, BUT IT */
/*     WILL LIKELY GENERATE AN IEEE FLOATING POINT UNDERFLOW ERROR */
/*     ON A SUN WORKSTATION.  REPLACE TENMAX WITH THE VALUE APPROPRIATE */
/*     FOR YOUR MACHINE. */

	tenmax = 320.;
	itnmax = 1;
	dnum = .1f;
L10:
	++itnmax;
	dnum *= .1;
	if (dnum > consts_1.zero) {
	    goto L10;
	}
	--itnmax;
	tenmax = (doublereal) itnmax;

	if (e1 > tenmax) {
	    x1 = tenmax;
	} else if (e1 < -tenmax) {
	    x1 = consts_1.zero;
	} else {
	    x1 = n1 * pow_dd(&consts_1.ten, &e1);
	}
	if (x2 != consts_1.zero) {
	    phi = atan2(x2, x1);
	} else {
	    phi = consts_1.zero;
	}
	d__1 = consts_1.half * (log(n3) + e3 * log(consts_1.ten));
	z__1.r = d__1, z__1.i = phi;
	c__->r = z__1.r, c__->i = z__1.i;
    }
    return 0;
} /* arydiv_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE EMULT                             * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Takes one base and exponent and multiplies it * */
/*     *    by another numbers base and exponent to give the product  * */
/*     *    in the form of base and exponent.                         * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************