pfqfn.cpp

Quantitative Finance C++ source code, asian option pricer

// testf2c2.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "f2c.h"


/* testpfq.f -- translated by f2c (version 20060506).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

//#include "f2c.h"


extern "C" integer i_dnnt(doublereal *);
extern "C" double d_imag(doublecomplex *);
extern "C" void z_div(doublecomplex *, doublecomplex *, doublecomplex *), pow_zz(
	    doublecomplex *, doublecomplex *, doublecomplex *);
extern "C" integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle(void);
extern "C" int s_stop(char *, ftnlen);
extern "C" double z_abs(doublecomplex *);
extern "C" integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

extern "C" void z_log(doublecomplex *, doublecomplex *);
extern "C" double exp(doublereal);
extern "C" void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
extern "C" double d_sign(doublereal *, doublereal *), pow_dd(doublereal *, 
	    doublereal *), z_abs(doublecomplex *);
extern "C"     double d_lg10(doublereal *), pow_di(doublereal *, integer *), d_int(
	    doublereal *), d_nint(doublereal *), d_imag(doublecomplex *);
extern "C"     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 *);
extern "C" double d_lg10(doublereal *), d_sign(doublereal *, doublereal *), pow_dd(
	    doublereal *, doublereal *), atan2(doublereal, doublereal), log(
	    doublereal);

extern "C"     double d_imag(doublecomplex *), z_abs(doublecomplex *), atan(doublereal), 
	    log(doublereal);
extern "C"     void z_log(doublecomplex *, doublecomplex *), z_div(doublecomplex *, 
	    doublecomplex *, doublecomplex *);
extern "C" double d_imag(doublecomplex *), atan(doublereal), pow_di(doublereal *, 
	    integer *), log(doublereal);
extern "C" integer i_dnnt(doublereal *), s_wsfe(cilist *), e_wsfe(void), do_fio(
	    integer *, char *, ftnlen);
    /* Subroutine */ int s_stop(char *, ftnlen);
extern "C" double sin(doublereal), sqrt(doublereal), atan2(doublereal, doublereal), 
	    exp(doublereal), cos(doublereal);


/* Common Block Declarations */

struct consts_1_ {
    doublereal zero, half, one, two, ten, eps;
};

#define consts_1 (*(struct consts_1_ *) &consts_)

struct {
    integer nout;
} io_;

#define io_1 io_

/* Initialized data */

struct {
    doublereal e_1[6];
    } consts_ = { 0., .5, 1., 2., 10., 1e-10 };


/* Table of constant values */

static integer c__9 = 9;
static integer c__1 = 1;
static integer c__25000 = 25000;
static doublereal c_b64 = 10.;

/*     **************************************************************** */
/*     *                                                              * */
/*     *    SOLUTION TO THE GENERALIZED HYPERGEOMETRIC FUNCTION       * */
/*     *                                                              * */
/*     *                           by                                 * */
/*     *                                                              * */
/*     *                      W. F. PERGER,                           * */
/*     *                                                              * */
/*     *              MARK NARDIN  and ATUL BHALLA                    * */
/*     *                                                              * */
/*     *                                                              * */
/*     *            Electrical Engineering Department                 * */
/*     *            Michigan Technological University                 * */
/*     *                  1400 Townsend Drive                         * */
/*     *                Houghton, MI  49931-1295   USA                * */
/*     *                     Copyright 1993                           * */
/*     *                                                              * */
/*     *               e-mail address: [e-mail]                    * */
/*     *                                                              * */
/*     *  Description : A numerical evaluator for the generalized     * */
/*     *    hypergeometric function for complex arguments with large  * */
/*     *    magnitudes using a direct summation of the Gauss series.  * */
/*     *    The method used allows an accuracy of up to thirteen      * */
/*     *    decimal places through the use of large integer arrays    * */
/*     *    and a single final division.                              * */
/*     *    (original subroutines for the confluent hypergeometric    * */
/*     *    written by Mark Nardin, 1989; modifications made to cal-  * */
/*     *    culate the generalized hypergeometric function were       * */
/*     *    written by W.F. Perger and A. Bhalla, June, 1990)         * */
/*     *                                                              * */
/*     *  The evaluation of the pFq series is accomplished by a func- * */
/*     *  ion call to PFQ, which is a double precision complex func-  * */
/*     *  tion.  The required input is:                               * */
/*     *  1. Double precision complex arrays A and B.  These are the  * */
/*     *     arrays containing the parameters in the numerator and de-* */
/*     *     nominator, respectively.                                 * */
/*     *  2. Integers IP and IQ.  These integers indicate the number  * */
/*     *     of numerator and denominator terms, respectively (these  * */
/*     *     are p and q in the pFq function).                        * */
/*     *  3. Double precision complex argument Z.                     * */
/*     *  4. Integer LNPFQ.  This integer should be set to '1' if the * */
/*     *     result from PFQ is to be returned as the natural logaritm* */
/*     *     of the series, or '0' if not.  The user can generally set* */
/*     *     LNPFQ = '0' and change it if required.                   * */
/*     *  5. Integer IX.  This integer should be set to '0' if the    * */
/*     *     user desires the program PFQ to estimate the number of   * */
/*     *     array terms (in A and B) to be used, or an integer       * */
/*     *     greater than zero specifying the number of integer pos-  * */
/*     *     itions to be used.  This input parameter is escpecially  * */
/*     *     useful as a means to check the results of a given run.   * */
/*     *     Specificially, if the user obtains a result for a given  * */
/*     *     set of parameters, then changes IX and re-runs the eval- * */
/*     *     uator, and if the number of array positions was insuffi- * */
/*     *     cient, then the two results will likely differ.  The rec-* */
/*     *     commended would be to generally set IX = '0' and then set* */
/*     *     it to 100 or so for a second run.  Note that the LENGTH  * */
/*     *     parameter currently sets the upper limit on IX to 777,   * */
/*     *     but that can easily be changed (it is a single PARAMETER * */
/*     *     statement) and the program recompiled.                   * */
/*     *  6. Integer NSIGFIG.  This integer specifies the requested   * */
/*     *     number of significant figures in the final result.  If   * */
/*     *     the user attempts to request more than the number of bits* */
/*     *     in the mantissa allows, the program will abort with an   * */
/*     *     appropriate error message.  The recommended value is 10. * */
/*     *                                                              * */
/*     *     Note: The variable NOUT is the file to which error mess- * */
/*     *           ages are written (default is 6).  This can be      * */
/*     *           changed in the FUNCTION PFQ to accomodate re-      * */
/*     *           of output to another file                          * */
/*     *                                                              * */
/*     *  Subprograms called: HYPER.                                  * */
/*     *                                                              * */
/*     **************************************************************** */

/* Double Complex */ VOID pfq_(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,\002IP=\002,1i2,3x,\002IQ=\002,1i2,3x,"
	    "\002AND ABS(Z)=\002,1e12.5,2x,/,\002 WHICH IS GREATER THAN ONE--"
	    "SERIES DOES\002,\002 NOT CONVERGE\002)";

    /* System generated locals */
    doublereal d__1;
    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6, z__7, z__8, z__9, z__10,
	     z__11, z__12, z__13, z__14, z__15, z__16, z__17;

    /* Builtin functions */
    //integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    //e_wsle(void);
    /* Subroutine */ int s_stop(char *, ftnlen);
   // double z_abs(doublecomplex *);
    //integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
    //void z_div(doublecomplex *, doublecomplex *, doublecomplex *), pow_zz(
	  //  doublecomplex *, doublecomplex *, doublecomplex *);
    //double d_imag(doublecomplex *);
    //integer i_dnnt(doublereal *);

    /* Local variables */
    static integer i__;
    static doublecomplex a1[2], b1[1], z1, gam1, gam2, gam3, gam4, gam5, gam6,
	     gam7;
    static doublereal diff, argi, argr, dnum;
    extern /* Double Complex */ VOID hyper_(doublecomplex *, doublecomplex *, 
	    doublecomplex *, integer *, integer *, doublecomplex *, integer *,
	     integer *, integer *);
    static doublecomplex hyper1, hyper2;
    extern /* Double Complex */ VOID cgamma_(doublecomplex *, doublecomplex *,
	     integer *);
    static doublereal precis;

    /* Fortran I/O blocks */
    static cilist io___1 = { 0, 0, 0, 0, 0 };
    static cilist io___2 = { 0, 0, 0, fmt_300, 0 };


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

    /* Function Body */
    io_1.nout = 6;
    z__1.r = consts_1.zero, z__1.i = consts_1.zero;
    if (z__->r == z__1.r && z__->i == z__1.i) {
	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 ;
    }
    if (*lnpfq != 0 && *lnpfq != 1) {
	io___1.ciunit = io_1.nout;
	s_wsle(&io___1);
	do_lio(&c__9, &c__1, " ERROR IN INPUT ARGUMENTS: LNPFQ /= 0 OR 1", (
		ftnlen)42);
	e_wsle();
	s_stop("", (ftnlen)0);
    }
    if (*ip > *iq && z_abs(z__) > consts_1.one && *ip > 2) {
	io___2.ciunit = io_1.nout;
	s_wsfe(&io___2);
	do_fio(&c__1, (char *)&(*ip), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&(*iq), (ftnlen)sizeof(integer));
	d__1 = z_abs(z__);
	do_fio(&c__1, (char *)&d__1, (ftnlen)sizeof(doublereal));
	e_wsfe();
	s_stop("", (ftnlen)0);
    }

/* For the 2F1 function with |z| > 1, use Abramowitz and Stegun, 15.3.7. */

    if (*ip == 2 && *iq == 1 && z_abs(z__) > consts_1.one) {
	z__2.r = consts_1.one, z__2.i = 0.;
	z_div(&z__1, &z__2, z__);
	z1.r = z__1.r, z1.i = z__1.i;
	a1[0].r = a[1].r, a1[0].i = a[1].i;
	z__2.r = consts_1.one - b[1].r, z__2.i = -b[1].i;
	z__1.r = z__2.r + a[1].r, z__1.i = z__2.i + a[1].i;
	a1[1].r = z__1.r, a1[1].i = z__1.i;
	z__2.r = consts_1.one - a[2].r, z__2.i = -a[2].i;
	z__1.r = z__2.r + a[1].r, z__1.i = z__2.i + a[1].i;
	b1[0].r = z__1.r, b1[0].i = z__1.i;
	cgamma_(&z__1, &b[1], lnpfq);
	gam1.r = z__1.r, gam1.i = z__1.i;
	z__2.r = a[2].r - a[1].r, z__2.i = a[2].i - a[1].i;
	cgamma_(&z__1, &z__2, lnpfq);
	gam2.r = z__1.r, gam2.i = z__1.i;
	cgamma_(&z__1, &a[2], lnpfq);
	gam3.r = z__1.r, gam3.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);
	gam4.r = z__1.r, gam4.i = z__1.i;
	z__2.r = a[1].r - a[2].r, z__2.i = a[1].i - a[2].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;
	z__2.r = b[1].r - a[2].r, z__2.i = b[1].i - a[2].i;
	cgamma_(&z__1, &z__2, lnpfq);
	gam7.r = z__1.r, gam7.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;
	a1[0].r = a[2].r, a1[0].i = a[2].i;
	z__2.r = consts_1.one - b[1].r, z__2.i = -b[1].i;
	z__1.r = z__2.r + a[2].r, z__1.i = z__2.i + a[2].i;
	a1[1].r = z__1.r, a1[1].i = z__1.i;
	z__2.r = consts_1.one - a[1].r, z__2.i = -a[1].i;
	z__1.r = z__2.r + a[2].r, z__1.i = z__2.i + a[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__5.r = gam1.r * gam2.r - gam1.i * gam2.i, z__5.i = gam1.r * gam2.i 
		+ gam1.i * gam2.r;
	z__7.r = -z__->r, z__7.i = -z__->i;
	z__8.r = -a[1].r, z__8.i = -a[1].i;
	pow_zz(&z__6, &z__7, &z__8);
	z__4.r = z__5.r * z__6.r - z__5.i * z__6.i, z__4.i = z__5.r * z__6.i 
		+ z__5.i * z__6.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__9.r = gam3.r * gam4.r - gam3.i * gam4.i, z__9.i = gam3.r * gam4.i 
		+ gam3.i * gam4.r;
	z_div(&z__2, &z__3, &z__9);
	z__13.r = gam1.r * gam5.r - gam1.i * gam5.i, z__13.i = gam1.r * 
		gam5.i + gam1.i * gam5.r;
	z__15.r = -z__->r, z__15.i = -z__->i;
	z__16.r = -a[2].r, z__16.i = -a[2].i;
	pow_zz(&z__14, &z__15, &z__16);
	z__12.r = z__13.r * z__14.r - z__13.i * z__14.i, z__12.i = z__13.r * 
		z__14.i + z__13.i * z__14.r;
	z__11.r = z__12.r * hyper2.r - z__12.i * hyper2.i, z__11.i = z__12.r *
		 hyper2.i + z__12.i * hyper2.r;
	z__17.r = gam6.r * gam7.r - gam6.i * gam7.i, z__17.i = gam6.r * 
		gam7.i + gam6.i * gam7.r;
	z_div(&z__10, &z__11, &z__17);
	z__1.r = z__2.r + z__10.r, z__1.i = z__2.i + z__10.i;
	 ret_val->r = z__1.r,  ret_val->i = z__1.i;
	return ;
    }
    if (*ip == 2 && *iq == 1 && z_abs(z__) > .9f) {
	if (*lnpfq == 1) {
	    goto L30;
	}

/*     Check to see if the Gamma function arguments are o.k.; if not, */
/*     then the series will have to be used. */

/*     PRECIS - MACHINE PRECISION */

	precis = consts_1.one;
L10:
	precis /= consts_1.two;
	dnum = precis + consts_1.one;
	if (dnum > consts_1.one) {
	    goto L10;
	}
	precis = consts_1.two * precis;
	for (i__ = 1; i__ <= 6; ++i__) {
	    if (i__ == 1) {
		argi = d_imag(&b[1]);
		argr = b[1].r;
	    } else if (i__ == 2) {
		z__2.r = b[1].r - a[1].r, z__2.i = b[1].i - a[1].i;
		z__1.r = z__2.r - a[2].r, z__1.i = z__2.i - a[2].i;
		argi = d_imag(&z__1);
		z__2.r = b[1].r - a[1].r, z__2.i = b[1].i - a[1].i;
		z__1.r = z__2.r - a[2].r, z__1.i = z__2.i - a[2].i;
		argr = z__1.r;
	    } else if (i__ == 3) {
		z__1.r = b[1].r - a[1].r, z__1.i = b[1].i - a[1].i;
		argi = d_imag(&z__1);
		z__1.r = b[1].r - a[1].r, z__1.i = b[1].i - a[1].i;
		argr = z__1.r;
	    } else if (i__ == 4) {
		z__2.r = a[1].r + a[2].r, z__2.i = a[1].i + a[2].i;
		z__1.r = z__2.r - b[1].r, z__1.i = z__2.i - b[1].i;
		argi = d_imag(&z__1);
		z__2.r = a[1].r + a[2].r, z__2.i = a[1].i + a[2].i;
		z__1.r = z__2.r - b[1].r, z__1.i = z__2.i - b[1].i;
		argr = z__1.r;
	    } else if (i__ == 5) {
		argi = d_imag(&a[1]);
		argr = a[1].r;
	    } else if (i__ == 6) {
		argi = d_imag(&a[2]);
		argr = a[2].r;
	    }

/*     CASES WHERE THE ARGUMENT IS REAL */

	    if (argi == consts_1.zero) {

/*     CASES WHERE THE ARGUMENT IS REAL AND NEGATIVE */

		if (argr <= consts_1.zero) {

/*     USE THE SERIES EXPANSION IF THE ARGUMENT IS TOO NEAR A POLE */

		    diff = (d__1 = (doublereal) i_dnnt(&argr) - argr, abs(
			    d__1));