pfqfn.cpp

Quantitative Finance C++ source code, asian option pricer

		    if (diff <= consts_1.two * precis) {
			goto L30;
		    }
		}
	    }
/* L20: */
	}
	cgamma_(&z__1, &b[1], lnpfq);
	gam1.r = z__1.r, gam1.i = z__1.i;
	z__3.r = b[1].r - a[1].r, z__3.i = b[1].i - a[1].i;
	z__2.r = z__3.r - a[2].r, z__2.i = z__3.i - a[2].i;
	cgamma_(&z__1, &z__2, lnpfq);
	gam2.r = z__1.r, gam2.i = z__1.i;
	z__2.r = b[1].r - a[1].r, z__2.i = b[1].i - a[1].i;
	cgamma_(&z__1, &z__2, lnpfq);
	gam3.r = z__1.r, gam3.i = z__1.i;
	z__2.r = b[1].r - a[2].r, z__2.i = b[1].i - a[2].i;
	cgamma_(&z__1, &z__2, lnpfq);
	gam4.r = z__1.r, gam4.i = z__1.i;
	z__3.r = a[1].r + a[2].r, z__3.i = a[1].i + a[2].i;
	z__2.r = z__3.r - b[1].r, z__2.i = z__3.i - b[1].i;
	cgamma_(&z__1, &z__2, lnpfq);
	gam5.r = z__1.r, gam5.i = z__1.i;
	cgamma_(&z__1, &a[1], lnpfq);
	gam6.r = z__1.r, gam6.i = z__1.i;
	cgamma_(&z__1, &a[2], lnpfq);
	gam7.r = z__1.r, gam7.i = z__1.i;
	a1[0].r = a[1].r, a1[0].i = a[1].i;
	a1[1].r = a[2].r, a1[1].i = a[2].i;
	z__3.r = a[1].r + a[2].r, z__3.i = a[1].i + a[2].i;
	z__2.r = z__3.r - b[1].r, z__2.i = z__3.i - b[1].i;
	z__1.r = z__2.r + consts_1.one, z__1.i = z__2.i;
	b1[0].r = z__1.r, b1[0].i = z__1.i;
	z__1.r = consts_1.one - z__->r, z__1.i = -z__->i;
	z1.r = z__1.r, z1.i = z__1.i;
	hyper_(&z__1, a1, b1, ip, iq, &z1, lnpfq, ix, nsigfig);
	hyper1.r = z__1.r, hyper1.i = z__1.i;
	z__1.r = b[1].r - a[1].r, z__1.i = b[1].i - a[1].i;
	a1[0].r = z__1.r, a1[0].i = z__1.i;
	z__1.r = b[1].r - a[2].r, z__1.i = b[1].i - a[2].i;
	a1[1].r = z__1.r, a1[1].i = z__1.i;
	z__3.r = b[1].r - a[1].r, z__3.i = b[1].i - a[1].i;
	z__2.r = z__3.r - a[2].r, z__2.i = z__3.i - a[2].i;
	z__1.r = z__2.r + consts_1.one, z__1.i = z__2.i;
	b1[0].r = z__1.r, b1[0].i = z__1.i;
	hyper_(&z__1, a1, b1, ip, iq, &z1, lnpfq, ix, nsigfig);
	hyper2.r = z__1.r, hyper2.i = z__1.i;
	z__4.r = gam1.r * gam2.r - gam1.i * gam2.i, z__4.i = gam1.r * gam2.i 
		+ gam1.i * gam2.r;
	z__3.r = z__4.r * hyper1.r - z__4.i * hyper1.i, z__3.i = z__4.r * 
		hyper1.i + z__4.i * hyper1.r;
	z__5.r = gam3.r * gam4.r - gam3.i * gam4.i, z__5.i = gam3.r * gam4.i 
		+ gam3.i * gam4.r;
	z_div(&z__2, &z__3, &z__5);
	z__11.r = consts_1.one - z__->r, z__11.i = -z__->i;
	z__13.r = b[1].r - a[1].r, z__13.i = b[1].i - a[1].i;
	z__12.r = z__13.r - a[2].r, z__12.i = z__13.i - a[2].i;
	pow_zz(&z__10, &z__11, &z__12);
	z__9.r = z__10.r * gam1.r - z__10.i * gam1.i, z__9.i = z__10.r * 
		gam1.i + z__10.i * gam1.r;
	z__8.r = z__9.r * gam5.r - z__9.i * gam5.i, z__8.i = z__9.r * gam5.i 
		+ z__9.i * gam5.r;
	z__7.r = z__8.r * hyper2.r - z__8.i * hyper2.i, z__7.i = z__8.r * 
		hyper2.i + z__8.i * hyper2.r;
	z__14.r = gam6.r * gam7.r - gam6.i * gam7.i, z__14.i = gam6.r * 
		gam7.i + gam6.i * gam7.r;
	z_div(&z__6, &z__7, &z__14);
	z__1.r = z__2.r + z__6.r, z__1.i = z__2.i + z__6.i;
	 ret_val->r = z__1.r,  ret_val->i = z__1.i;
	return ;
    }
L30:
    hyper_(&z__1, &a[1], &b[1], ip, iq, z__, lnpfq, ix, nsigfig);
     ret_val->r = z__1.r,  ret_val->i = z__1.i;
    return ;
} /* pfq_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                   FUNCTION BITS                              * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Determines the number of significant figures  * */
/*     *    of machine precision to arrive at the size of the array   * */
/*     *    the numbers must be stored in to get the accuracy of the  * */
/*     *    solution.                                                 * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************************** */
doublereal bits_(void)
{
    /* System generated locals */
    doublereal ret_val;

    /* Local variables */
    static doublereal bit, bit2;
    static integer count;

    bit = 1.f;
    count = 0;
L10:
    ++count;
    bit2 = bit * 2.f;
    bit = bit2 + 1.f;
    if (bit - bit2 != 0.f) {
	goto L10;
    }
    ret_val = (doublereal) (count - 3);
    return ret_val;
} /* bits_ */

/*     **************************************************************** */
/*     *                                                              * */
/*     *                   FUNCTION HYPER                             * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Function that sums the Gauss series.          * */
/*     *                                                              * */
/*     *  Subprograms called: ARMULT, ARYDIV, BITS, CMPADD, CMPMUL,   * */
/*     *                      IPREMAX.                                * */
/*     *                                                              * */
/*     **************************************************************** */
/* Double Complex */ VOID hyper_(doublecomplex * ret_val, doublecomplex *a, 
	doublecomplex *b, integer *ip, integer *iq, doublecomplex *z__, 
	integer *lnpfq, integer *ix, integer *nsigfig)
{
    /* Format strings */
    static char fmt_300[] = "(1x,\002WARNING - REAL PART OF A(\002,1i2,\002)"
	    " WAS SET TO ZERO\002)";
    static char fmt_301[] = "(1x,\002WARNING - IMAG PART OF A(\002,1i2,\002)"
	    " WAS SET TO ZERO\002)";
    static char fmt_302[] = "(1x,\002WARNING - REAL PART OF B(\002,1i2,\002)"
	    " WAS SET TO ZERO\002)";
    static char fmt_303[] = "(1x,\002WARNING - IMAG PART OF B(\002,1i2,\002)"
	    " WAS SET TO ZERO\002)";
    static char fmt_304[] = "(1x,\002ERROR - ARGUMENT B(\002,1i2,\002) WAS E"
	    "QUAL TO ZERO\002)";
    static char fmt_305[] = "(1x,\002ERROR - ARGUMENT B(\002,1i2,\002) WAS A"
	    " NEGATIVE\002,\002 INTEGER\002)";
    static char fmt_306[] = "(1x,\002ERROR IN FN HYPER: L MUST BE < \002,1i4)"
	    ;
    static char fmt_307[] = "(1x,\002 WARNING--THE NUMBER OF SIGNIFICANT FIG"
	    "URES REQU\002,\002ESTED\002,/,\002IS GREATER THAN THE MACHINE PR"
	    "ECISION--\002,\002FINAL ANSWER\002,/,\002WILL BE ACCURATE TO ONLY"
	    "\002,i3,\002 DIGITS\002)";

    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5;
    doublereal d__1, d__2, d__3;
    doublecomplex z__1, z__2, z__3, z__4, z__5;

    /* Builtin functions */
//    double d_lg10(doublereal *), pow_di(doublereal *, integer *), d_int(
//	    doublereal *), d_nint(doublereal *), d_imag(doublecomplex *);
//    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
//	     s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
//	    e_wsle(void), i_dnnt(doublereal *);
    ///* Subroutine */ int s_stop(char *, ftnlen);
//    void z_log(doublecomplex *, doublecomplex *);
//    double exp(doublereal);
//    void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
//    double d_sign(doublereal *, doublereal *), pow_dd(doublereal *, 
//	    doublereal *), z_abs(doublecomplex *);

    /* Local variables */
    static integer i__, l;
    static doublereal x;
    static integer i1;
    static doublereal ai[10], ci[10], ar[10], cr[10], xi, xl, xr, ai2[10], 
	    ci2[10], ar2[10], cr2[10], qi1[25002], qi2[25002], wk1[25002], 
	    qr1[25002], qr2[25002], wk2[25002], wk3[25002], wk4[25002], wk5[
	    25002], wk6[25002], mx1, mx2, xr2, xi2;
    static integer ii10, ir10;
    static doublereal ri10, cnt, rr10, log2, dum1, dum2, accy;
    static integer ibit, lmax;
    extern doublereal bits_(void);
    static doublecomplex temp;
    static doublereal rmax, numi[25002], sumi[25002];
    static integer rexp;
    static doublereal numr[25002], sumr[25002];
    static doublecomplex cdum1, cdum2, temp1;
    static integer nmach;
    static doublereal creal;
    static doublecomplex final;
    static integer ixcnt;
    static doublereal expon;
    extern /* Subroutine */ int cmpadd_(doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, integer *, doublereal *);
    static doublereal sigfig, denomi[25002];
    extern /* Double Complex */ VOID factor_(doublecomplex *, doublecomplex *)
	    ;
    static doublereal denomr[25002];
    extern /* Subroutine */ int cmpmul_(doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, doublereal *, integer *, doublereal *), arydiv_(
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublecomplex *, integer *, integer *, doublereal *, integer *);
    static integer icount;
    extern /* Subroutine */ int armult_(doublereal *, doublereal *, 
	    doublereal *, doublereal *, integer *, doublereal *);
    static doublecomplex oldtemp;
    extern integer ipremax_(doublecomplex *, doublecomplex *, integer *, 
	    integer *, doublecomplex *);

    /* Fortran I/O blocks */
    static cilist io___41 = { 0, 0, 0, fmt_300, 0 };
    static cilist io___42 = { 0, 0, 0, fmt_301, 0 };
    static cilist io___43 = { 0, 0, 0, fmt_302, 0 };
    static cilist io___44 = { 0, 0, 0, fmt_303, 0 };
    static cilist io___45 = { 0, 0, 0, 0, 0 };
    static cilist io___46 = { 0, 0, 0, 0, 0 };
    static cilist io___49 = { 0, 0, 0, fmt_304, 0 };
    static cilist io___50 = { 0, 0, 0, fmt_305, 0 };
    static cilist io___59 = { 0, 0, 0, fmt_306, 0 };
    static cilist io___60 = { 0, 0, 0, fmt_307, 0 };


    /* Parameter adjustments */
    --a;
    --b;

    /* Function Body */
    consts_1.zero = 0.;
    log2 = d_lg10(&consts_1.two);
    ibit = (integer) bits_();
    i__1 = ibit / 2;
    rmax = pow_di(&consts_1.two, &i__1);
    i__1 = ibit / 4;
    sigfig = pow_di(&consts_1.two, &i__1);

    i__1 = *ip;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = i1;
	ar2[i1 - 1] = a[i__2].r * sigfig;
	ar[i1 - 1] = d_int(&ar2[i1 - 1]);
	d__1 = (ar2[i1 - 1] - ar[i1 - 1]) * rmax;
	ar2[i1 - 1] = d_nint(&d__1);
	ai2[i1 - 1] = d_imag(&a[i1]) * sigfig;
	ai[i1 - 1] = d_int(&ai2[i1 - 1]);
	d__1 = (ai2[i1 - 1] - ai[i1 - 1]) * rmax;
	ai2[i1 - 1] = d_nint(&d__1);
/* L10: */
    }
    i__1 = *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = i1;
	cr2[i1 - 1] = b[i__2].r * sigfig;
	cr[i1 - 1] = d_int(&cr2[i1 - 1]);
	d__1 = (cr2[i1 - 1] - cr[i1 - 1]) * rmax;
	cr2[i1 - 1] = d_nint(&d__1);
	ci2[i1 - 1] = d_imag(&b[i1]) * sigfig;
	ci[i1 - 1] = d_int(&ci2[i1 - 1]);
	d__1 = (ci2[i1 - 1] - ci[i1 - 1]) * rmax;
	ci2[i1 - 1] = d_nint(&d__1);
/* L20: */
    }
    xr2 = z__->r * sigfig;
    xr = d_int(&xr2);
    d__1 = (xr2 - xr) * rmax;
    xr2 = d_nint(&d__1);
    xi2 = d_imag(z__) * sigfig;
    xi = d_int(&xi2);
    d__1 = (xi2 - xi) * rmax;
    xi2 = d_nint(&d__1);

/*     WARN THE USER THAT THE INPUT VALUE WAS SO CLOSE TO ZERO THAT IT */
/*     WAS SET EQUAL TO ZERO. */

    i__1 = *ip;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = i1;
	if (a[i__2].r != consts_1.zero && ar[i1 - 1] == consts_1.zero && ar2[
		i1 - 1] == consts_1.zero) {
	    io___41.ciunit = io_1.nout;
	    s_wsfe(&io___41);
	    do_fio(&c__1, (char *)&i1, (ftnlen)sizeof(integer));
	    e_wsfe();
	}
	if (d_imag(&a[i1]) != consts_1.zero && ai[i1 - 1] == consts_1.zero && 
		ai2[i1 - 1] == consts_1.zero) {
	    io___42.ciunit = io_1.nout;
	    s_wsfe(&io___42);
	    do_fio(&c__1, (char *)&i1, (ftnlen)sizeof(integer));
	    e_wsfe();
	}
/* L30: */
    }
    i__1 = *iq;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = i1;
	if (b[i__2].r != consts_1.zero && cr[i1 - 1] == consts_1.zero && cr2[
		i1 - 1] == consts_1.zero) {
	    io___43.ciunit = io_1.nout;
	    s_wsfe(&io___43);
	    do_fio(&c__1, (char *)&i1, (ftnlen)sizeof(integer));
	    e_wsfe();
	}
	if (d_imag(&b[i1]) != consts_1.zero && ci[i1 - 1] == consts_1.zero && 
		ci2[i1 - 1] == consts_1.zero) {
	    io___44.ciunit = io_1.nout;
	    s_wsfe(&io___44);
	    do_fio(&c__1, (char *)&i1, (ftnlen)sizeof(integer));
	    e_wsfe();
	}
/* L40: */
    }
    if (z__->r != consts_1.zero && xr == consts_1.zero && xr2 == 
	    consts_1.zero) {
	io___45.ciunit = io_1.nout;
	s_wsle(&io___45);
	do_lio(&c__9, &c__1, " WARNING - REAL PART OF Z WAS SET TO ZERO", (
		ftnlen)41);
	e_wsle();
	d__1 = d_imag(z__);
	z__1.r = consts_1.zero, z__1.i = d__1;
	z__->r = z__1.r, z__->i = z__1.i;
    }
    if (d_imag(z__) != consts_1.zero && xi == consts_1.zero && xi2 == 
	    consts_1.zero) {
	io___46.ciunit = io_1.nout;
	s_wsle(&io___46);
	do_lio(&c__9, &c__1, " WARNING - IMAG PART OF Z WAS SET TO ZERO", (
		ftnlen)41);
	e_wsle();
	d__1 = z__->r;
	z__1.r = d__1, z__1.i = consts_1.zero;
	z__->r = z__1.r, z__->i = z__1.i;
    }

/*     SCREENING OF NUMERATOR ARGUMENTS FOR NEGATIVE INTEGERS OR ZERO. */
/*     ICOUNT WILL FORCE THE SERIES TO TERMINATE CORRECTLY. */

    i__1 = (integer) bits_();
    d__1 = pow_di(&consts_1.two, &i__1);
    nmach = (integer) d_lg10(&d__1);
    icount = -1;
    i__1 = *ip;
    for (i1 = 1; i1 <= i__1; ++i1) {
	if (ar2[i1 - 1] == consts_1.zero && ar[i1 - 1] == consts_1.zero && 
		ai2[i1 - 1] == consts_1.zero && ai[i1 - 1] == consts_1.zero) {
	    z__1.r = consts_1.one, z__1.i = consts_1.zero;
	     ret_val->r = z__1.r,  ret_val->i = z__1.i;
	    return ;
	}
	i__2 = i1;
	if (ai[i1 - 1] == consts_1.zero && ai2[i1 - 1] == consts_1.zero && a[
		i__2].r < consts_1.zero) {
	    i__2 = i1;
	    i__3 = i1;
	    d__2 = a[i__3].r;
	    i__4 = -nmach;
	    if ((d__1 = a[i__2].r - (doublereal) i_dnnt(&d__2), abs(d__1)) < 
		    pow_di(&consts_1.ten, &i__4)) {
		if (icount != -1) {
/* Computing MIN */
		    i__4 = i1;
		    d__1 = a[i__4].r;
		    i__2 = icount, i__3 = -i_dnnt(&d__1);
		    icount = min(i__2,i__3);
		} else {
		    i__2 = i1;
		    d__1 = a[i__2].r;
		    icount = -i_dnnt(&d__1);