TOMS_707.C

Go to the documentation of this file.

 
/* toms_707.f -- translated by f2c (version 19980516).
   You must link the resulting object file with the libraries:
        -lf2c -lm   (in that order)
*/
 
/* $Id: TOMS_707.C,v 1.2.2.11 2019/05/28 11:13:02 garloff Exp $ */
 
#if defined(__GNUC__) && __GNUC__ == 2 && defined(__cplusplus)
# undef __cplusplus     /* workaround for buggy gcc-2.96 */
#endif
 
#include "tbci/specfun/prototypes.h"
#include "tbci/specfun/prototypes2.h"
 
#ifdef __cplusplus
extern "C" {
#endif
#include "tbci/lapack/f2c.h"
 
/* Common Block Declarations */
 
struct {
    doublereal y;
} stcom_;
 
#define stcom_1 stcom_
 
/* Table of constant values */
 
static doublereal c_b8 = 2.;
static doublereal c_b53 = 1.;
static doublereal c_b65 = 10.;
static integer c__7 = 7;
static integer c__1 = 1;
static integer c__9 = 9;
static integer c__3 = 3;
static integer c__2 = 2;
 
/*     **************************************************************** */
int aradd_(const doublereal *a, const doublereal *b,
        doublereal *c__, const integer *l, const doublereal *rmax);
int arsub_(const doublereal *a, const doublereal *b,
        doublereal *c__, const integer *l, const doublereal *rmax);
 
int armult_(const doublereal *a, const doublereal *b, doublereal *c__,
        const integer *l, const doublereal *rmax);
int arydiv_(const doublereal *ar, const doublereal *ai,
        const doublereal *br, const doublereal *bi,
        doublecomplex *c__, const integer *l, const integer *lnchf,
        const doublereal *rmax, const integer *bit);
 
int cmpadd_(const doublereal *ar, const doublereal *ai,
        const doublereal *br, const doublereal *bi,
        doublereal *cr, doublereal *ci,
        const integer *l, const doublereal *rmax);
int cmpmul_(const doublereal *ar, const doublereal *ai,
        const doublereal *br, const doublereal *bi,
        doublereal *cr, doublereal *ci,
        const integer *l, const doublereal *rmax);
 
    extern integer bits_(void);
    extern doublereal store_(doublereal *x);
 
    extern /* Subroutine */ int conv12_(doublecomplex *cn, doublereal *cae), conv21_(doublereal *cae, doublecomplex *cn);
 
    extern /* Subroutine */ int emult_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
    extern /* Subroutine */ int ecpmul_(doublereal *a, doublereal *b, doublereal *c__), eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
    extern /* Subroutine */ int ecpdiv_(doublereal *a, doublereal *b, doublereal *c__);
    extern /* Subroutine */ int ediv_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
    extern /* Subroutine */ int eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
    extern /* Subroutine */ int esub_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
 
 
 
/* Double Complex */
VOID chgf_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z__, const integer *l, const integer *lnchf);
 
/*      ALGORITHM 707, COLLECTED ALGORITHMS FROM ACM. */
/*      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, */
/*      VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 345-349. */
/*     **************************************************************** */
/*     *                                                              * */
/*     *      SOLUTION TO THE CONFLUENT HYPERGEOMETRIC FUNCTION       * */
/*     *                                                              * */
/*     *                           by                                 * */
/*     *                                                              * */
/*     *                      MARK NARDIN,                            * */
/*     *                                                              * */
/*     *              W. F. PERGER and ATUL BHALLA                    * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Michigan Technological University, Copyright 1989           * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : A numerical evaluator for the confluent       * */
/*     *    hypergeometric function for complex arguments with large  * */
/*     *    magnitudes using a direct summation of the Kummer series. * */
/*     *    The method used allows an accuracy of up to thirteen      * */
/*     *    decimal places through the use of large real arrays       * */
/*     *    and a single final division.  LNCHF is a variable which   * */
/*     *    selects how the result should be represented.  A '0' will * */
/*     *    return the value in standard exponential form.  A '1'     * */
/*     *    will return the LOG of the result.  IP is an integer      * */
/*     *    variable that specifies how many array positions are      * */
/*     *    desired (usually 10 is sufficient).  Setting IP=0 causes  * */
/*     *    the program to estimate the number of array positions.    * */
/*     *                                                              * */
/*     *    The confluent hypergeometric function is the solution to  * */
/*     *    the differential equation:                                * */
/*     *                                                              * */
/*     *             zf"(z) + (a-z)f'(z) - bf(z) = 0                  * */
/*     *                                                              * */
/*     *  Subprograms called: BITS, CHGF                              * */
/*     *                                                              * */
/*     **************************************************************** */
/* Double Complex */ /*VOID*/
void conhyp_(doublecomplex *ret_val,
        const doublecomplex *a, const doublecomplex *b,
        const doublecomplex *z__, const integer *lnchf, const integer *ip)
{
    /* System generated locals */
    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6, z__7;
 
    /* Builtin functions: prototypes.h */
 
    /* Local variables */
    static doublereal term1, term2;
    static integer i__;
    static doublereal nterm, fx, ang, max__;
 
    if (z_abs(z__) != 0.) {
        ang = atan2(d_imag(z__), (doublereal) z__->r);
    } else {
        ang = 1.;
    }
    if (abs(ang) < 1.570796) {
        ang = 1.;
    } else {
        ang = sin(abs(ang) - 1.570796325) + 1.;
    }
    max__ = 0.;
    nterm = 0.;
    fx = 0.;
    term1 = 0.;
L10:
    nterm += 1;
    z__4.r = a->r + nterm, z__4.i = a->i;
    z__3.r = z__4.r - 1, z__3.i = z__4.i;
    z__2.r = z__3.r * z__->r - z__3.i * z__->i, z__2.i = z__3.r * z__->i +
            z__3.i * z__->r;
    z__7.r = b->r + nterm, z__7.i = b->i;
    z__6.r = z__7.r - 1, z__6.i = z__7.i;
    z__5.r = nterm * z__6.r, z__5.i = nterm * z__6.i;
    z_div(&z__1, &z__2, &z__5);
    term2 = z_abs(&z__1);
    if (term2 == 0.) {
        goto L20;
    }
    if (term2 < 1.) {
        if (a->r + nterm - 1 > 1.) {
            if (b->r + nterm - 1 > 1.) {
                if (term2 - term1 < 0.) {
                    goto L20;
                }
            }
        }
    }
    fx += log(term2);
    if (fx > max__) {
        max__ = fx;
    }
    term1 = term2;
    goto L10;
L20:
    max__ = max__ * 2 / (bits_() * .69314718056);
    i__ = (integer) (max__ * ang) + 7;
    if (i__ < 5) {
        i__ = 5;
    }
    if (*ip > i__) {
        i__ = *ip;
    }
    chgf_(&z__1, a, b, z__, &i__, lnchf);
     ret_val->r = z__1.r,  ret_val->i = z__1.i;
    return ;
} /* conhyp_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                   FUNCTION BITS                              * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Determines the number of significant figures  * */
/*     *    of machine precision to arrive at the size of the array   * */
/*     *    the numbers must must be stored in to get the accuracy    * */
/*     *    of the solution.                                          * */
/*     *                                                              * */
/*     *  Subprogram called: STORE                                    * */
/*     *                                                              * */
/*     **************************************************************** */
integer bits_(void)
{
    /* System generated locals */
    integer ret_val;
    doublereal d__1;
 
    /* Local variables */
    static integer count;
    static doublereal bit, bit2;
 
    bit = (float)1.;
    count = 0;
L10:
    ++count;
    d__1 = bit * (float)2.;
    bit2 = store_(&d__1);
    d__1 = bit2 + (float)1.;
    bit = store_(&d__1);
    if (bit - bit2 != (float)0.) {
        goto L10;
    }
    ret_val = count;
    return ret_val;
} /* bits_ */
 
doublereal store_(doublereal *x)
{
    /* System generated locals */
    doublereal ret_val;
 
 
/* *********************************************************** */
 
 
/*   This function forces its argument X to be stored in a */
/* memory location, thus providing a means of determining */
/* floating point number characteristics (such as the machine */
/* precision) when it is necessary to avoid computation in */
/* high precision registers. */
 
/* On input: */
 
/*       X = Value to be stored. */
 
/* X is not altered by this function. */
 
/* On output: */
 
/*       STORE = Value of X after it has been stored and */
/*               possibly truncated or rounded to the double */
/*               precision word length. */
 
/* Modules required by STORE:  None */
 
/* *********************************************************** */
 
    stcom_1.y = *x;
    ret_val = stcom_1.y;
    return ret_val;
} /* store_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                   FUNCTION CHGF                              * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Function that sums the Kummer series and      * */
/*     *    returns the solution of the confluent hypergeometric      * */
/*     *    function.                                                 * */
/*     *                                                              * */
/*     *  Subprograms called: ARMULT, ARYDIV, BITS, CMPADD, CMPMUL    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Double Complex */ VOID chgf_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z__, const integer *l, const integer *lnchf)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1, d__2;
    doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;
 
    /* Builtin functions */
 
    /* Local variables */
    static doublereal rmax, numi[779], sumi[779], numr[779], sumr[779];
    static integer i__;
    static doublecomplex final;
    static doublereal ai, ci, ar, cr;
    static doublereal xi, sigfig, xr, denomi[779], denomr[779], ai2;
    static doublereal ci2;
    static doublereal ar2, cr2, qi1[779], qi2[779], xi2, qr1[779], qr2[779],
            mx1, mx2, xr2;
    static integer bit;
    static doublereal cnt;
 
    bit = bits_();
    i__1 = bit / 2;
    rmax = pow_di(&c_b8, &i__1);
    i__1 = bit / 4;
    sigfig = pow_di(&c_b8, &i__1);
    ar2 = a->r * sigfig;
    ar = d_int(&ar2);
    d__1 = (ar2 - ar) * rmax;
    ar2 = d_nint(&d__1);
    ai2 = d_imag(a) * sigfig;
    ai = d_int(&ai2);
    d__1 = (ai2 - ai) * rmax;
    ai2 = d_nint(&d__1);
    cr2 = b->r * sigfig;
    cr = d_int(&cr2);
    d__1 = (cr2 - cr) * rmax;
    cr2 = d_nint(&d__1);
    ci2 = d_imag(b) * sigfig;
    ci = d_int(&ci2);
    d__1 = (ci2 - ci) * rmax;
    ci2 = d_nint(&d__1);
    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);
    sumr[0] = 1.;
    sumi[0] = 1.;
    numr[0] = 1.;
    numi[0] = 1.;
    denomr[0] = 1.;
    denomi[0] = 1.;
    i__1 = *l + 1;
    for (i__ = 0; i__ <= i__1; ++i__) {
        sumr[i__ + 1] = 0.;
        sumi[i__ + 1] = 0.;
        numr[i__ + 1] = 0.;
        numi[i__ + 1] = 0.;
        denomr[i__ + 1] = 0.;
        denomi[i__ + 1] = 0.;
/* L100: */
    }
    sumr[2] = 1.;
    numr[2] = 1.;
    denomr[2] = 1.;
    cnt = sigfig;
L110:
    if (sumr[2] < (float).5) {
        mx1 = sumi[*l + 2];
    } else if (sumi[2] < (float).5) {
        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] < (float).5) {
        mx2 = numi[*l + 2];
    } else if (numi[2] < (float).5) {
        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 > (float)2.) {
        if (cr > 0.) {
            z__3.r = ar, z__3.i = ai;
            z__4.r = xr, z__4.i = xi;
            z__2.r = z__3.r * z__4.r - z__3.i * z__4.i, z__2.i = z__3.r *
                    z__4.i + z__3.i * z__4.r;
            z__6.r = cr, z__6.i = ci;
            z__5.r = cnt * z__6.r, z__5.i = cnt * z__6.i;
            z_div(&z__1, &z__2, &z__5);
            if (z_abs(&z__1) <= 1.) {
                goto L190;
            }
        }
    }
    cmpmul_(sumr, sumi, &cr, &ci, qr1, qi1, l, &rmax);
    cmpmul_(sumr, sumi, &cr2, &ci2, qr2, qi2, l, &rmax);
    qr2[*l + 2] += -1;
    qi2[*l + 2] += -1;
    cmpadd_(qr1, qi1, qr2, qi2, sumr, sumi, l, &rmax);
    armult_(sumr, &cnt, sumr, l, &rmax);
    armult_(sumi, &cnt, sumi, l, &rmax);
    cmpmul_(denomr, denomi, &cr, &ci, qr1, qi1, l, &rmax);
    cmpmul_(denomr, denomi, &cr2, &ci2, qr2, qi2, l, &rmax);
    qr2[*l + 2] += -1;
    qi2[*l + 2] += -1;
    cmpadd_(qr1, qi1, qr2, qi2, denomr, denomi, l, &rmax);
    armult_(denomr, &cnt, denomr, l, &rmax);
    armult_(denomi, &cnt, denomi, l, &rmax);
    cmpmul_(numr, numi, &ar, &ai, qr1, qi1, l, &rmax);
    cmpmul_(numr, numi, &ar2, &ai2, qr2, qi2, l, &rmax);
    qr2[*l + 2] += -1;
    qi2[*l + 2] += -1;
    cmpadd_(qr1, qi1, qr2, qi2, numr, numi, l, &rmax);
    cmpmul_(numr, numi, &xr, &xi, qr1, qi1, l, &rmax);
    cmpmul_(numr, numi, &xr2, &xi2, qr2, qi2, l, &rmax);
    qr2[*l + 2] += -1;
    qi2[*l + 2] += -1;
    cmpadd_(qr1, qi1, qr2, qi2, numr, numi, l, &rmax);
    cmpadd_(sumr, sumi, numr, numi, sumr, sumi, l, &rmax);
    cnt += sigfig;
    ar += sigfig;
    cr += sigfig;
    goto L110;
L190:
    arydiv_(sumr, sumi, denomr, denomi, &final, l, lnchf, &rmax, &bit);
     ret_val->r = final.r,  ret_val->i = final.i;
    return ;
} /* chgf_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 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_(const doublereal *a, const doublereal *b,
        doublereal *c__, const integer *l, const doublereal *rmax)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1;
 
    /* Builtin functions */
 
    /* Local variables */
    static integer i__, j, ediff;
    static doublereal z__[779];
 
    /* Parameter adjustments */
    ++c__;
    ++b;
    ++a;
 
    /* Function Body */
    i__1 = *l + 1;
    for (i__ = 0; i__ <= i__1; ++i__) {
        z__[i__ + 1] = 0.;
/* L110: */
    }
    d__1 = a[*l + 1] - b[*l + 1];
    ediff = (integer) d_nint(&d__1);
    if (abs(a[1]) < (float).5 || ediff <= -(*l)) {
        goto L111;
    }
    if (abs(b[1]) < (float).5 || ediff >= *l) {
        goto L113;
    }
    goto L115;
L111:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
        c__[i__] = b[i__];
/* L112: */
    }
    goto L311;
L113:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
        c__[i__] = a[i__];
/* L114: */
    }
    goto L311;
L115:
    z__[0] = a[-1];
    if ((d__1 = a[-1] - b[-1], abs(d__1)) < (float).5) {
        goto L200;
    }
    if (ediff > 0) {
        z__[*l + 2] = a[*l + 1];
        goto L233;
    }
    if (ediff < 0) {
        z__[*l + 2] = b[*l + 1];
        z__[0] = b[-1];
        goto L266;
    }
    i__1 = *l;
    for (i__ = 1; i__ <= i__1; ++i__) {
        if (a[i__] > b[i__]) {
            z__[*l + 2] = a[*l + 1];
            goto L233;
        }
        if (a[i__] < b[i__]) {
            z__[*l + 2] = b[*l + 1];
            z__[0] = b[-1];
            goto L266;
        }
/* L120: */
    }
    goto L300;
L200:
    if (ediff > 0) {
        goto L203;
    }
    if (ediff < 0) {
        goto L207;
    }
    z__[*l + 2] = a[*l + 1];
    for (i__ = *l; i__ >= 1; --i__) {
        z__[i__ + 1] = a[i__] + b[i__] + z__[i__ + 1];
        if (z__[i__ + 1] >= *rmax) {
            z__[i__ + 1] -= *rmax;
            z__[i__] = 1.;
        }
/* L201: */
    }
    if (z__[1] > (float).5) {
        for (i__ = *l; i__ >= 1; --i__) {
            z__[i__ + 1] = z__[i__];
/* L202: */
        }
        z__[*l + 2] += 1.;
        z__[1] = 0.;
    }
    goto L300;
L203:
    z__[*l + 2] = a[*l + 1];
    i__1 = ediff + 1;
    for (i__ = *l; i__ >= i__1; --i__) {
        z__[i__ + 1] = a[i__] + b[i__ - ediff] + z__[i__ + 1];
        if (z__[i__ + 1] >= *rmax) {
            z__[i__ + 1] -= *rmax;
            z__[i__] = 1.;
        }
/* L204: */
    }
    for (i__ = ediff; i__ >= 1; --i__) {
        z__[i__ + 1] = a[i__] + z__[i__ + 1];
        if (z__[i__ + 1] >= *rmax) {
            z__[i__ + 1] -= *rmax;
            z__[i__] = 1.;
        }
/* L205: */
    }
    if (z__[1] > (float).5) {
        for (i__ = *l; i__ >= 1; --i__) {
            z__[i__ + 1] = z__[i__];
/* L206: */
        }
        z__[*l + 2] += 1;
        z__[1] = 0.;
    }
    goto L300;
L207:
    z__[*l + 2] = b[*l + 1];
    i__1 = 1 - ediff;
    for (i__ = *l; i__ >= i__1; --i__) {
        z__[i__ + 1] = a[i__ + ediff] + b[i__] + z__[i__ + 1];
        if (z__[i__ + 1] >= *rmax) {
            z__[i__ + 1] -= *rmax;
            z__[i__] = 1.;
        }
/* L208: */
    }
    for (i__ = -ediff; i__ >= 1; --i__) {
        z__[i__ + 1] = b[i__] + z__[i__ + 1];
        if (z__[i__ + 1] >= *rmax) {
            z__[i__ + 1] -= *rmax;
            z__[i__] = 1.;
        }
/* L209: */
    }
    if (z__[1] > (float).5) {
        for (i__ = *l; i__ >= 1; --i__) {
            z__[i__ + 1] = z__[i__];
/* L210: */
        }
        z__[*l + 2] += 1.;
        z__[1] = 0.;
    }
    goto L300;
L233:
    if (ediff > 0) {
        goto L243;
    }
    for (i__ = *l; i__ >= 1; --i__) {
        z__[i__ + 1] = a[i__] - b[i__] + z__[i__ + 1];
        if (z__[i__ + 1] < 0.) {
            z__[i__ + 1] += *rmax;
            z__[i__] = -1.;
        }
/* L234: */
    }
    goto L290;
L243:
    i__1 = ediff + 1;
    for (i__ = *l; i__ >= i__1; --i__) {
        z__[i__ + 1] = a[i__] - b[i__ - ediff] + z__[i__ + 1];
        if (z__[i__ + 1] < 0.) {
            z__[i__ + 1] += *rmax;
            z__[i__] = -1.;
        }
/* L244: */
    }
    for (i__ = ediff; i__ >= 1; --i__) {
        z__[i__ + 1] = a[i__] + z__[i__ + 1];
        if (z__[i__ + 1] < 0.) {
            z__[i__ + 1] += *rmax;
            z__[i__] = -1.;
        }
/* L245: */
    }
    goto L290;
L266:
    if (ediff < 0) {
        goto L276;
    }
    for (i__ = *l; i__ >= 1; --i__) {
        z__[i__ + 1] = b[i__] - a[i__] + z__[i__ + 1];
        if (z__[i__ + 1] < 0.) {
            z__[i__ + 1] += *rmax;
            z__[i__] = -1.;
        }
/* L267: */
    }
    goto L290;
L276:
    i__1 = 1 - ediff;
    for (i__ = *l; i__ >= i__1; --i__) {
        z__[i__ + 1] = b[i__] - a[i__ + ediff] + z__[i__ + 1];
        if (z__[i__ + 1] < 0.) {
            z__[i__ + 1] += *rmax;
            z__[i__] = -1.;
        }
/* L277: */
    }
    for (i__ = -ediff; i__ >= 1; --i__) {
        z__[i__ + 1] = b[i__] + z__[i__ + 1];
        if (z__[i__ + 1] < 0.) {
            z__[i__ + 1] += *rmax;
            z__[i__] = -1.;
        }
/* L278: */
    }
L290:
    if (z__[2] > (float).5) {
        goto L300;
    }
    i__ = 1;
L291:
    ++i__;
    if (z__[i__ + 1] < (float).5 && i__ < *l + 1) {
        goto L291;
    }
    if (i__ == *l + 1) {
        z__[0] = 1.;
        z__[*l + 2] = 0.;
        goto L300;
    }
/* L292: */
    i__1 = *l + 1 - i__;
    for (j = 1; j <= i__1; ++j) {
        z__[j + 1] = z__[j + i__];
/* L293: */
    }
    i__1 = *l;
    for (j = *l + 2 - i__; j <= i__1; ++j) {
        z__[j + 1] = 0.;
/* L294: */
    }
    z__[*l + 2] = z__[*l + 2] - i__ + 1;
L300:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
        c__[i__] = z__[i__ + 1];
/* L310: */
    }
L311:
    if (c__[1] < (float).5) {
        c__[-1] = 1.;
        c__[*l + 1] = 0.;
    }
    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_(const doublereal *a, const doublereal *b,
        doublereal *c__, const integer *l, const doublereal *rmax)
{
    /* System generated locals */
    integer i__1;
 
    /* Local variables */
    static integer i__;
    static doublereal b2[779];
 
    /* Parameter adjustments */
    ++c__;
    ++b;
    ++a;
 
    /* Function Body */
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
        b2[i__ + 1] = b[i__];
/* L100: */
    }
    b2[0] *= -1.;
    aradd_(&a[-1], b2, &c__[-1], l, rmax);
    return 0;
} /* arsub_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ARMULT                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  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_(const doublereal *a, const doublereal *b, doublereal *c__, const integer *l, const doublereal *rmax)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1;
 
    /* Builtin functions */
 
    /* Local variables */
    static doublereal rmax2;
    static integer i__;
    static doublereal z__[779], carry, b2;
 
    /* Parameter adjustments */
    ++c__;
    ++a;
 
    /* Function Body */
    rmax2 = 1. / *rmax;
    z__[0] = d_sign(&c_b53, b) * a[-1];
    b2 = abs(*b);
    z__[*l + 2] = a[*l + 1];
    i__1 = *l;
    for (i__ = 0; i__ <= i__1; ++i__) {
        z__[i__ + 1] = 0.;
/* L100: */
    }
    if (b2 <= 1e-10 || a[1] <= 1e-10) {
        z__[0] = 1.;
        z__[*l + 2] = 0.;
        goto L198;
    }
    for (i__ = *l; i__ >= 1; --i__) {
        z__[i__ + 1] = a[i__] * b2 + z__[i__ + 1];
        if (z__[i__ + 1] >= *rmax) {
            d__1 = z__[i__ + 1] / *rmax;
            carry = d_int(&d__1);
            z__[i__ + 1] -= carry * *rmax;
            z__[i__] = carry;
        }
/* L110: */
    }
    if (z__[1] < (float).5) {
        goto L150;
    }
    for (i__ = *l; i__ >= 1; --i__) {
        z__[i__ + 1] = z__[i__];
/* L120: */
    }
    z__[*l + 2] += 1.;
    z__[1] = 0.;
L150:
L198:
    i__1 = *l + 1;
    for (i__ = -1; i__ <= i__1; ++i__) {
        c__[i__] = z__[i__ + 1];
/* L199: */
    }
    if (c__[1] < (float).5) {
        c__[-1] = 1.;
        c__[*l + 1] = 0.;
    }
    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_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublereal *cr, doublereal *ci, const integer *l, const doublereal *rmax)
{
 
    /* Parameter adjustments */
    ++ci;
    ++cr;
    ++bi;
    ++br;
    ++ai;
    ++ar;
 
    /* Function Body */
    aradd_(&ar[-1], &br[-1], &cr[-1], l, rmax);
    aradd_(&ai[-1], &bi[-1], &ci[-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, integer *l, doublereal *rmax)
{
    /* Parameter adjustments */
    ++ci;
    ++cr;
    ++bi;
    ++br;
    ++ai;
    ++ar;
 
    /* Function Body */
    arsub_(&ar[-1], &br[-1], &cr[-1], l, rmax);
    arsub_(&ai[-1], &bi[-1], &ci[-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_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublereal *cr, doublereal *ci, const integer *l, const doublereal *rmax)
{
    static doublereal d1[779], d2[779];
 
    /* Parameter adjustments */
    ++ci;
    ++cr;
    ++ai;
    ++ar;
 
    /* Function Body */
    armult_(&ar[-1], br, d1, l, rmax);
    armult_(&ai[-1], bi, d2, l, rmax);
    arsub_(d1, d2, &cr[-1], l, rmax);
    armult_(&ar[-1], bi, d1, l, rmax);
    armult_(&ai[-1], br, d2, l, rmax);
    aradd_(d1, d2, &ci[-1], l, rmax);
    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  * */
/*     *    two different forms.  Either standard scientific or as    * */
/*     *    the log of the number.                                    * */
/*     *                                                              * */
/*     *  Subprograms called: CONV21, CONV12, EADD, ECPDIV, EMULT     * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int arydiv_(const doublereal *ar, const doublereal *ai, const doublereal *br, const doublereal *bi, doublecomplex *c__, const integer *l, const integer *lnchf, const doublereal *rmax, const integer *bit)
{
    /* System generated locals */
    doublereal d__1;
    doublecomplex z__1;
 
    /* Builtin functions */
 
    /* Local variables */
    static integer rexp;
    static doublereal x;
    static doublereal e1, e2, e3;
    static doublereal 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;
    /* Parameter adjustments */
    ++bi;
    ++br;
    ++ai;
    ++ar;
 
    /* Function Body */
    rexp = *bit / 2;
    x = rexp * (ar[*l + 1] - 2);
    rr10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
    ir10 = (integer) rr10;
    rr10 -= ir10;
    x = rexp * (ai[*l + 1] - 2);
    ri10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
    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_b65, &rr10);
    dum2 *= pow_dd(&c_b65, &ri10);
    z__1.r = dum1, z__1.i = dum2;
    conv12_(&z__1, ae);
    ae[2] += ir10;
    ae[3] += ii10;
    x = rexp * (br[*l + 1] - 2);
    rr10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
    ir10 = (integer) rr10;
    rr10 -= ir10;
    x = rexp * (bi[*l + 1] - 2);
    ri10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
    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_b65, &rr10);
    dum2 *= pow_dd(&c_b65, &ri10);
    z__1.r = dum1, z__1.i = dum2;
    conv12_(&z__1, be);
    be[2] += ir10;
    be[3] += ii10;
    ecpdiv_(ae, be, ce);
    if (*lnchf == 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];
        if (e1 > 74.) {
            x1 = 1e75;
        } else if (e1 < -74.) {
            x1 = 0.;
        } else {
            x1 = n1 * pow_dd(&c_b65, &e1);
        }
        phi = atan2(x2, x1);
        d__1 = (log(n3) + e3 * log(10.)) * .5;
        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                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int emult_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
{
    *nf = *n1 * *n2;
    *ef = *e1 + *e2;
    if (abs(*nf) >= 10.) {
        *nf /= 10.;
        *ef += 1.;
    }
    return 0;
} /* emult_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE EDIV                              * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : returns the solution in the form of base and  * */
/*     *    exponent of the division of two exponential numbers.      * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int ediv_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
{
    *nf = *n1 / *n2;
    *ef = *e1 - *e2;
    if (abs(*nf) < 1. && *nf != 0.) {
        *nf *= 10.;
        *ef += -1.;
    }
    return 0;
} /* ediv_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE EADD                              * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Returns the sum of two numbers in the form    * */
/*     *    of a base and an exponent.                                * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
{
    /* Builtin functions */
 
    /* Local variables */
    static doublereal ediff;
 
    ediff = *e1 - *e2;
    if (ediff > 36.) {
        *nf = *n1;
        *ef = *e1;
    } else if (ediff < -36.) {
        *nf = *n2;
        *ef = *e2;
    } else {
        *nf = *n1 * pow_dd(&c_b65, &ediff) + *n2;
        *ef = *e2;
L400:
        if (abs(*nf) < 10.) {
            goto L410;
        }
        *nf /= 10.;
        *ef += 1.;
        goto L400;
L410:
        if (abs(*nf) >= 1. || *nf == 0.) {
            goto L420;
        }
        *nf *= 10.;
        *ef += -1.;
        goto L410;
    }
L420:
    return 0;
} /* eadd_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ESUB                              * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Returns the solution to the subtraction of    * */
/*     *    two numbers in the form of base and exponent.             * */
/*     *                                                              * */
/*     *  Subprograms called: EADD                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int esub_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
{
    /* System generated locals */
    doublereal d__1;
 
    /* Local variables */
 
    d__1 = *n2 * -1.;
    eadd_(n1, e1, &d__1, e2, nf, ef);
    return 0;
} /* esub_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE CONV12                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Converts a number from complex notation to a  * */
/*     *    form of a 2x2 real array.                                 * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int conv12_(doublecomplex *cn, doublereal *cae)
{
    /* Builtin functions */
 
    /* Parameter adjustments */
    cae -= 3;
 
    /* Function Body */
    cae[3] = cn->r;
    cae[5] = 0.;
L300:
    if (abs(cae[3]) < 10.) {
        goto L310;
    }
    cae[3] /= 10.;
    cae[5] += 1.;
    goto L300;
L310:
    if (abs(cae[3]) >= 1. || cae[3] == 0.) {
        goto L320;
    }
    cae[3] *= 10.;
    cae[5] += -1.;
    goto L310;
L320:
    cae[4] = d_imag(cn);
    cae[6] = 0.;
L330:
    if (abs(cae[4]) < 10.) {
        goto L340;
    }
    cae[4] /= 10.;
    cae[6] += 1.;
    goto L330;
L340:
    if (abs(cae[4]) >= 1. || cae[4] == 0.) {
        goto L350;
    }
    cae[4] *= 10.;
    cae[6] += -1.;
    goto L340;
L350:
    return 0;
} /* conv12_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE CONV21                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Converts a number represented in a 2x2 real   * */
/*     *    array to the form of a complex number.                    * */
/*     *                                                              * */
/*     *  Subprograms called: none                                    * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int conv21_(doublereal *cae, doublecomplex *cn)
{
    /* System generated locals */
    doublereal d__1, d__2;
    doublecomplex z__1;
 
    /* Builtin functions */
 
    /* Parameter adjustments */
    cae -= 3;
 
    /* Function Body */
    if (cae[5] > 75. || cae[6] > 75.) {
        cn->r = 1e75, cn->i = 1e75;
    } else if (cae[6] < -75.) {
        d__1 = cae[3] * pow_dd(&c_b65, &cae[5]);
        z__1.r = d__1, z__1.i = 0.;
        cn->r = z__1.r, cn->i = z__1.i;
    } else {
        d__1 = cae[3] * pow_dd(&c_b65, &cae[5]);
        d__2 = cae[4] * pow_dd(&c_b65, &cae[6]);
        z__1.r = d__1, z__1.i = d__2;
        cn->r = z__1.r, cn->i = z__1.i;
    }
    return 0;
} /* conv21_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ECPMUL                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Multiplies two numbers which are each         * */
/*     *    represented in the form of a two by two array and returns * */
/*     *    the solution in the same form.                            * */
/*     *                                                              * */
/*     *  Subprograms called: EMULT, ESUB, EADD                       * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int ecpmul_(doublereal *a, doublereal *b, doublereal *c__)
{
    static doublereal c2[4]     /* was [2][2] */, e1, e2;
    static doublereal n1, n2;
    /* Parameter adjustments */
    c__ -= 3;
    b -= 3;
    a -= 3;
 
    /* Function Body */
    emult_(&a[3], &a[5], &b[3], &b[5], &n1, &e1);
    emult_(&a[4], &a[6], &b[4], &b[6], &n2, &e2);
    esub_(&n1, &e1, &n2, &e2, c2, &c2[2]);
    emult_(&a[3], &a[5], &b[4], &b[6], &n1, &e1);
    emult_(&a[4], &a[6], &b[3], &b[5], &n2, &e2);
    eadd_(&n1, &e1, &n2, &e2, &c__[4], &c__[6]);
    c__[3] = c2[0];
    c__[5] = c2[2];
    return 0;
} /* ecpmul_ */
 
/*     **************************************************************** */
/*     *                                                              * */
/*     *                 SUBROUTINE ECPDIV                            * */
/*     *                                                              * */
/*     *                                                              * */
/*     *  Description : Divides two numbers and returns the solution. * */
/*     *    All numbers are represented by a 2x2 array.               * */
/*     *                                                              * */
/*     *  Subprograms called: EADD, ECPMUL, EDIV, EMULT               * */
/*     *                                                              * */
/*     **************************************************************** */
/* Subroutine */ int ecpdiv_(doublereal *a, doublereal *b, doublereal *c__)
{
    static doublereal b2[4]     /* was [2][2] */, c2[4] /* was [2][2] */, e1,
            e2, e3;
    static doublereal n1, n2, n3;
    /* Parameter adjustments */
    c__ -= 3;
    b -= 3;
    a -= 3;
 
    /* Function Body */
    b2[0] = b[3];
    b2[2] = b[5];
    b2[1] = b[4] * -1.;
    b2[3] = b[6];
    ecpmul_(&a[3], b2, c2);
    emult_(&b[3], &b[5], &b[3], &b[5], &n1, &e1);
    emult_(&b[4], &b[6], &b[4], &b[6], &n2, &e2);
    eadd_(&n1, &e1, &n2, &e2, &n3, &e3);
    ediv_(c2, &c2[2], &n3, &e3, &c__[3], &c__[5]);
    ediv_(&c2[1], &c2[3], &n3, &e3, &c__[4], &c__[6]);
    return 0;
} /* ecpdiv_ */
 
/* Main program */ int MAIN__(void)
{
    /* Format strings */
    static char fmt_300[] = "(1x,\002 RESULT FROM CONHYP=\002,1p,2d25.12)";
    static char fmt_301[] = "(1x,\002 EXPECTED RESULT=\002,3x,1p,2d25.12)";
 
    /* System generated locals */
    doublecomplex z__1;
    olist o__1;
 
    /* Builtin functions */
 
    /* Local variables */
    static doublecomplex a, b, z__;
    static integer lnchf;
    static doublereal ai, ar;
    static integer ip;
    static doublecomplex chf;
 
    /* Fortran I/O blocks */
    static cilist io___88 = { 0, 5, 0, 0, 0 };
    static cilist io___90 = { 0, 6, 0, 0, 0 };
    static cilist io___91 = { 0, 5, 0, 0, 0 };
    static cilist io___93 = { 0, 6, 0, 0, 0 };
    static cilist io___94 = { 0, 5, 0, 0, 0 };
    static cilist io___96 = { 0, 6, 0, 0, 0 };
    static cilist io___97 = { 0, 5, 0, 0, 0 };
    static cilist io___99 = { 0, 6, 0, 0, 0 };
    static cilist io___100 = { 0, 5, 0, 0, 0 };
    static cilist io___102 = { 0, 6, 0, 0, 0 };
    static cilist io___104 = { 0, 6, 0, fmt_300, 0 };
    static cilist io___107 = { 0, 6, 0, fmt_301, 0 };
 
 
 
    o__1.oerr = 0;
    o__1.ounit = 5;
    o__1.ofnmlen = 4;
    o__1.ofnm = "DATA";
    o__1.orl = 0;
    o__1.osta = "OLD";
    o__1.oacc = 0;
    o__1.ofm = 0;
    o__1.oblnk = 0;
    f_open(&o__1);
    s_rsle(&io___88);
    do_lio(&c__7, &c__1, (char *)&a, (ftnlen)sizeof(doublecomplex));
    e_rsle();
    s_wsle(&io___90);
    do_lio(&c__9, &c__1, " A= ", (ftnlen)4);
    do_lio(&c__7, &c__1, (char *)&a, (ftnlen)sizeof(doublecomplex));
    e_wsle();
    s_rsle(&io___91);
    do_lio(&c__7, &c__1, (char *)&b, (ftnlen)sizeof(doublecomplex));
    e_rsle();
    s_wsle(&io___93);
    do_lio(&c__9, &c__1, " B= ", (ftnlen)4);
    do_lio(&c__7, &c__1, (char *)&b, (ftnlen)sizeof(doublecomplex));
    e_wsle();
    s_rsle(&io___94);
    do_lio(&c__7, &c__1, (char *)&z__, (ftnlen)sizeof(doublecomplex));
    e_rsle();
    s_wsle(&io___96);
    do_lio(&c__9, &c__1, " Z= ", (ftnlen)4);
    do_lio(&c__7, &c__1, (char *)&z__, (ftnlen)sizeof(doublecomplex));
    e_wsle();
    s_rsle(&io___97);
    do_lio(&c__3, &c__1, (char *)&lnchf, (ftnlen)sizeof(integer));
    e_rsle();
    s_wsle(&io___99);
    do_lio(&c__9, &c__1, " LNCHF= ", (ftnlen)8);
    do_lio(&c__3, &c__1, (char *)&lnchf, (ftnlen)sizeof(integer));
    e_wsle();
    s_rsle(&io___100);
    do_lio(&c__3, &c__1, (char *)&ip, (ftnlen)sizeof(integer));
    e_rsle();
    s_wsle(&io___102);
    do_lio(&c__9, &c__1, " IP= ", (ftnlen)5);
    do_lio(&c__3, &c__1, (char *)&ip, (ftnlen)sizeof(integer));
    e_wsle();
    conhyp_(&z__1, &a, &b, &z__, &lnchf, &ip);
    chf.r = z__1.r, chf.i = z__1.i;
    s_wsfe(&io___104);
    do_fio(&c__2, (char *)&chf, (ftnlen)sizeof(doublereal));
    e_wsfe();
    ar = 2.31145634403113e-12;
    ai = -1.96169649634905e-11;
    s_wsfe(&io___107);
    do_fio(&c__1, (char *)&ar, (ftnlen)sizeof(doublereal));
    do_fio(&c__1, (char *)&ai, (ftnlen)sizeof(doublereal));
    e_wsfe();
    s_stop("", (ftnlen)0);
    return 0;
} /* MAIN__ */
 
/* Main program alias */ int sample_ (void) { return MAIN__ (); }
#ifdef __cplusplus
}
#endif