TBCI Numerical high perf. C++ Library: TOMS

00001 
00005 /* toms_707.f -- translated by f2c (version 19980516).
00006    You must link the resulting object file with the libraries:
00007         -lf2c -lm   (in that order)
00008 */
00009 
00010 /* $Id: TOMS_707.C,v 1.2.2.9 2008/07/17 17:07:08 garloff Exp $ */
00011 
00012 #if defined(__GNUC__) && __GNUC__ == 2 && defined(__cplusplus)
00013 # undef __cplusplus     /* workaround for buggy gcc-2.96 */
00014 #endif
00015 
00016 #include "prototypes.h"
00017 #include "prototypes2.h"
00018 
00019 #ifdef __cplusplus
00020 extern "C" {
00021 #endif
00022 #include "f2c.h"
00023 
00024 /* Common Block Declarations */
00025 
00026 struct {
00027     doublereal y;
00028 } stcom_;
00029 
00030 #define stcom_1 stcom_
00031 
00032 /* Table of constant values */
00033 
00034 static doublereal c_b8 = 2.;
00035 static doublereal c_b53 = 1.;
00036 static doublereal c_b65 = 10.;
00037 static integer c__7 = 7;
00038 static integer c__1 = 1;
00039 static integer c__9 = 9;
00040 static integer c__3 = 3;
00041 static integer c__2 = 2;
00042 
00043 /*     **************************************************************** */
00044 int aradd_(const doublereal *a, const doublereal *b, 
00045         doublereal *c__, const integer *l, const doublereal *rmax);
00046 int arsub_(const doublereal *a, const doublereal *b, 
00047         doublereal *c__, const integer *l, const doublereal *rmax);
00048 
00049 int armult_(const doublereal *a, const doublereal *b, doublereal *c__, 
00050         const integer *l, const doublereal *rmax);
00051 int arydiv_(const doublereal *ar, const doublereal *ai, 
00052         const doublereal *br, const doublereal *bi, 
00053         doublecomplex *c__, const integer *l, const integer *lnchf, 
00054         const doublereal *rmax, const integer *bit); 
00055 
00056 int cmpadd_(const doublereal *ar, const doublereal *ai, 
00057         const doublereal *br, const doublereal *bi, 
00058         doublereal *cr, doublereal *ci, 
00059         const integer *l, const doublereal *rmax);
00060 int cmpmul_(const doublereal *ar, const doublereal *ai, 
00061         const doublereal *br, const doublereal *bi, 
00062         doublereal *cr, doublereal *ci, 
00063         const integer *l, const doublereal *rmax);
00064 
00065     extern integer bits_(void);
00066     extern doublereal store_(doublereal *x);
00067     
00068     extern /* Subroutine */ int conv12_(doublecomplex *cn, doublereal *cae), conv21_(doublereal *cae, doublecomplex *cn);
00069 
00070     extern /* Subroutine */ int emult_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
00071     extern /* Subroutine */ int ecpmul_(doublereal *a, doublereal *b, doublereal *c__), eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
00072     extern /* Subroutine */ int ecpdiv_(doublereal *a, doublereal *b, doublereal *c__);
00073     extern /* Subroutine */ int ediv_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
00074     extern /* Subroutine */ int eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
00075     extern /* Subroutine */ int esub_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef);
00076 
00077 
00078 
00079 /* Double Complex */ 
00080 VOID chgf_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z__, const integer *l, const integer *lnchf);
00081 
00082 /*      ALGORITHM 707, COLLECTED ALGORITHMS FROM ACM. */
00083 /*      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, */
00084 /*      VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 345-349. */
00085 /*     **************************************************************** */
00086 /*     *                                                              * */
00087 /*     *      SOLUTION TO THE CONFLUENT HYPERGEOMETRIC FUNCTION       * */
00088 /*     *                                                              * */
00089 /*     *                           by                                 * */
00090 /*     *                                                              * */
00091 /*     *                      MARK NARDIN,                            * */
00092 /*     *                                                              * */
00093 /*     *              W. F. PERGER and ATUL BHALLA                    * */
00094 /*     *                                                              * */
00095 /*     *                                                              * */
00096 /*     *  Michigan Technological University, Copyright 1989           * */
00097 /*     *                                                              * */
00098 /*     *                                                              * */
00099 /*     *  Description : A numerical evaluator for the confluent       * */
00100 /*     *    hypergeometric function for complex arguments with large  * */
00101 /*     *    magnitudes using a direct summation of the Kummer series. * */
00102 /*     *    The method used allows an accuracy of up to thirteen      * */
00103 /*     *    decimal places through the use of large real arrays       * */
00104 /*     *    and a single final division.  LNCHF is a variable which   * */
00105 /*     *    selects how the result should be represented.  A '0' will * */
00106 /*     *    return the value in standard exponential form.  A '1'     * */
00107 /*     *    will return the LOG of the result.  IP is an integer      * */
00108 /*     *    variable that specifies how many array positions are      * */
00109 /*     *    desired (usually 10 is sufficient).  Setting IP=0 causes  * */
00110 /*     *    the program to estimate the number of array positions.    * */
00111 /*     *                                                              * */
00112 /*     *    The confluent hypergeometric function is the solution to  * */
00113 /*     *    the differential equation:                                * */
00114 /*     *                                                              * */
00115 /*     *             zf"(z) + (a-z)f'(z) - bf(z) = 0                  * */
00116 /*     *                                                              * */
00117 /*     *  Subprograms called: BITS, CHGF                              * */
00118 /*     *                                                              * */
00119 /*     **************************************************************** */
00120 /* Double Complex */ /*VOID*/
00121 void conhyp_(doublecomplex *ret_val, 
00122         const doublecomplex *a, const doublecomplex *b, 
00123         const doublecomplex *z__, const integer *lnchf, const integer *ip)
00124 {
00125     /* System generated locals */
00126     doublecomplex z__1, z__2, z__3, z__4, z__5, z__6, z__7;
00127 
00128     /* Builtin functions: prototypes.h */
00129 
00130     /* Local variables */
00131     static doublereal term1, term2;
00132     static integer i__;
00133     static doublereal nterm, fx, ang, max__;
00134 
00135     if (z_abs(z__) != 0.) {
00136         ang = atan2(d_imag(z__), (doublereal) z__->r);
00137     } else {
00138         ang = 1.;
00139     }
00140     if (abs(ang) < 1.570796) {
00141         ang = 1.;
00142     } else {
00143         ang = sin(abs(ang) - 1.570796325) + 1.;
00144     }
00145     max__ = 0.;
00146     nterm = 0.;
00147     fx = 0.;
00148     term1 = 0.;
00149 L10:
00150     nterm += 1;
00151     z__4.r = a->r + nterm, z__4.i = a->i;
00152     z__3.r = z__4.r - 1, z__3.i = z__4.i;
00153     z__2.r = z__3.r * z__->r - z__3.i * z__->i, z__2.i = z__3.r * z__->i + 
00154             z__3.i * z__->r;
00155     z__7.r = b->r + nterm, z__7.i = b->i;
00156     z__6.r = z__7.r - 1, z__6.i = z__7.i;
00157     z__5.r = nterm * z__6.r, z__5.i = nterm * z__6.i;
00158     z_div(&z__1, &z__2, &z__5);
00159     term2 = z_abs(&z__1);
00160     if (term2 == 0.) {
00161         goto L20;
00162     }
00163     if (term2 < 1.) {
00164         if (a->r + nterm - 1 > 1.) {
00165             if (b->r + nterm - 1 > 1.) {
00166                 if (term2 - term1 < 0.) {
00167                     goto L20;
00168                 }
00169             }
00170         }
00171     }
00172     fx += log(term2);
00173     if (fx > max__) {
00174         max__ = fx;
00175     }
00176     term1 = term2;
00177     goto L10;
00178 L20:
00179     max__ = max__ * 2 / (bits_() * .69314718056);
00180     i__ = (integer) (max__ * ang) + 7;
00181     if (i__ < 5) {
00182         i__ = 5;
00183     }
00184     if (*ip > i__) {
00185         i__ = *ip;
00186     }
00187     chgf_(&z__1, a, b, z__, &i__, lnchf);
00188      ret_val->r = z__1.r,  ret_val->i = z__1.i;
00189     return ;
00190 } /* conhyp_ */
00191 
00192 /*     **************************************************************** */
00193 /*     *                                                              * */
00194 /*     *                   FUNCTION BITS                              * */
00195 /*     *                                                              * */
00196 /*     *                                                              * */
00197 /*     *  Description : Determines the number of significant figures  * */
00198 /*     *    of machine precision to arrive at the size of the array   * */
00199 /*     *    the numbers must must be stored in to get the accuracy    * */
00200 /*     *    of the solution.                                          * */
00201 /*     *                                                              * */
00202 /*     *  Subprogram called: STORE                                    * */
00203 /*     *                                                              * */
00204 /*     **************************************************************** */
00205 integer bits_(void)
00206 {
00207     /* System generated locals */
00208     integer ret_val;
00209     doublereal d__1;
00210 
00211     /* Local variables */
00212     static integer count;
00213     static doublereal bit, bit2;
00214 
00215     bit = (float)1.;
00216     count = 0;
00217 L10:
00218     ++count;
00219     d__1 = bit * (float)2.;
00220     bit2 = store_(&d__1);
00221     d__1 = bit2 + (float)1.;
00222     bit = store_(&d__1);
00223     if (bit - bit2 != (float)0.) {
00224         goto L10;
00225     }
00226     ret_val = count;
00227     return ret_val;
00228 } /* bits_ */
00229 
00230 doublereal store_(doublereal *x)
00231 {
00232     /* System generated locals */
00233     doublereal ret_val;
00234 
00235 
00236 /* *********************************************************** */
00237 
00238 
00239 /*   This function forces its argument X to be stored in a */
00240 /* memory location, thus providing a means of determining */
00241 /* floating point number characteristics (such as the machine */
00242 /* precision) when it is necessary to avoid computation in */
00243 /* high precision registers. */
00244 
00245 /* On input: */
00246 
00247 /*       X = Value to be stored. */
00248 
00249 /* X is not altered by this function. */
00250 
00251 /* On output: */
00252 
00253 /*       STORE = Value of X after it has been stored and */
00254 /*               possibly truncated or rounded to the double */
00255 /*               precision word length. */
00256 
00257 /* Modules required by STORE:  None */
00258 
00259 /* *********************************************************** */
00260 
00261     stcom_1.y = *x;
00262     ret_val = stcom_1.y;
00263     return ret_val;
00264 } /* store_ */
00265 
00266 /*     **************************************************************** */
00267 /*     *                                                              * */
00268 /*     *                   FUNCTION CHGF                              * */
00269 /*     *                                                              * */
00270 /*     *                                                              * */
00271 /*     *  Description : Function that sums the Kummer series and      * */
00272 /*     *    returns the solution of the confluent hypergeometric      * */
00273 /*     *    function.                                                 * */
00274 /*     *                                                              * */
00275 /*     *  Subprograms called: ARMULT, ARYDIV, BITS, CMPADD, CMPMUL    * */
00276 /*     *                                                              * */
00277 /*     **************************************************************** */
00278 /* Double Complex */ VOID chgf_(doublecomplex *ret_val, const doublecomplex *a, const doublecomplex *b, const doublecomplex *z__, const integer *l, const integer *lnchf)
00279 {
00280     /* System generated locals */
00281     integer i__1;
00282     doublereal d__1, d__2;
00283     doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;
00284 
00285     /* Builtin functions */
00286 
00287     /* Local variables */
00288     static doublereal rmax, numi[779], sumi[779], numr[779], sumr[779];
00289     static integer i__;
00290     static doublecomplex final;
00291     static doublereal ai, ci, ar, cr;
00292     static doublereal xi, sigfig, xr, denomi[779], denomr[779], ai2;
00293     static doublereal ci2;
00294     static doublereal ar2, cr2, qi1[779], qi2[779], xi2, qr1[779], qr2[779], 
00295             mx1, mx2, xr2;
00296     static integer bit;
00297     static doublereal cnt;
00298 
00299     bit = bits_();
00300     i__1 = bit / 2;
00301     rmax = pow_di(&c_b8, &i__1);
00302     i__1 = bit / 4;
00303     sigfig = pow_di(&c_b8, &i__1);
00304     ar2 = a->r * sigfig;
00305     ar = d_int(&ar2);
00306     d__1 = (ar2 - ar) * rmax;
00307     ar2 = d_nint(&d__1);
00308     ai2 = d_imag(a) * sigfig;
00309     ai = d_int(&ai2);
00310     d__1 = (ai2 - ai) * rmax;
00311     ai2 = d_nint(&d__1);
00312     cr2 = b->r * sigfig;
00313     cr = d_int(&cr2);
00314     d__1 = (cr2 - cr) * rmax;
00315     cr2 = d_nint(&d__1);
00316     ci2 = d_imag(b) * sigfig;
00317     ci = d_int(&ci2);
00318     d__1 = (ci2 - ci) * rmax;
00319     ci2 = d_nint(&d__1);
00320     xr2 = z__->r * sigfig;
00321     xr = d_int(&xr2);
00322     d__1 = (xr2 - xr) * rmax;
00323     xr2 = d_nint(&d__1);
00324     xi2 = d_imag(z__) * sigfig;
00325     xi = d_int(&xi2);
00326     d__1 = (xi2 - xi) * rmax;
00327     xi2 = d_nint(&d__1);
00328     sumr[0] = 1.;
00329     sumi[0] = 1.;
00330     numr[0] = 1.;
00331     numi[0] = 1.;
00332     denomr[0] = 1.;
00333     denomi[0] = 1.;
00334     i__1 = *l + 1;
00335     for (i__ = 0; i__ <= i__1; ++i__) {
00336         sumr[i__ + 1] = 0.;
00337         sumi[i__ + 1] = 0.;
00338         numr[i__ + 1] = 0.;
00339         numi[i__ + 1] = 0.;
00340         denomr[i__ + 1] = 0.;
00341         denomi[i__ + 1] = 0.;
00342 /* L100: */
00343     }
00344     sumr[2] = 1.;
00345     numr[2] = 1.;
00346     denomr[2] = 1.;
00347     cnt = sigfig;
00348 L110:
00349     if (sumr[2] < (float).5) {
00350         mx1 = sumi[*l + 2];
00351     } else if (sumi[2] < (float).5) {
00352         mx1 = sumr[*l + 2];
00353     } else {
00354 /* Computing MAX */
00355         d__1 = sumr[*l + 2], d__2 = sumi[*l + 2];
00356         mx1 = max(d__1,d__2);
00357     }
00358     if (numr[2] < (float).5) {
00359         mx2 = numi[*l + 2];
00360     } else if (numi[2] < (float).5) {
00361         mx2 = numr[*l + 2];
00362     } else {
00363 /* Computing MAX */
00364         d__1 = numr[*l + 2], d__2 = numi[*l + 2];
00365         mx2 = max(d__1,d__2);
00366     }
00367     if (mx1 - mx2 > (float)2.) {
00368         if (cr > 0.) {
00369             z__3.r = ar, z__3.i = ai;
00370             z__4.r = xr, z__4.i = xi;
00371             z__2.r = z__3.r * z__4.r - z__3.i * z__4.i, z__2.i = z__3.r * 
00372                     z__4.i + z__3.i * z__4.r;
00373             z__6.r = cr, z__6.i = ci;
00374             z__5.r = cnt * z__6.r, z__5.i = cnt * z__6.i;
00375             z_div(&z__1, &z__2, &z__5);
00376             if (z_abs(&z__1) <= 1.) {
00377                 goto L190;
00378             }
00379         }
00380     }
00381     cmpmul_(sumr, sumi, &cr, &ci, qr1, qi1, l, &rmax);
00382     cmpmul_(sumr, sumi, &cr2, &ci2, qr2, qi2, l, &rmax);
00383     qr2[*l + 2] += -1;
00384     qi2[*l + 2] += -1;
00385     cmpadd_(qr1, qi1, qr2, qi2, sumr, sumi, l, &rmax);
00386     armult_(sumr, &cnt, sumr, l, &rmax);
00387     armult_(sumi, &cnt, sumi, l, &rmax);
00388     cmpmul_(denomr, denomi, &cr, &ci, qr1, qi1, l, &rmax);
00389     cmpmul_(denomr, denomi, &cr2, &ci2, qr2, qi2, l, &rmax);
00390     qr2[*l + 2] += -1;
00391     qi2[*l + 2] += -1;
00392     cmpadd_(qr1, qi1, qr2, qi2, denomr, denomi, l, &rmax);
00393     armult_(denomr, &cnt, denomr, l, &rmax);
00394     armult_(denomi, &cnt, denomi, l, &rmax);
00395     cmpmul_(numr, numi, &ar, &ai, qr1, qi1, l, &rmax);
00396     cmpmul_(numr, numi, &ar2, &ai2, qr2, qi2, l, &rmax);
00397     qr2[*l + 2] += -1;
00398     qi2[*l + 2] += -1;
00399     cmpadd_(qr1, qi1, qr2, qi2, numr, numi, l, &rmax);
00400     cmpmul_(numr, numi, &xr, &xi, qr1, qi1, l, &rmax);
00401     cmpmul_(numr, numi, &xr2, &xi2, qr2, qi2, l, &rmax);
00402     qr2[*l + 2] += -1;
00403     qi2[*l + 2] += -1;
00404     cmpadd_(qr1, qi1, qr2, qi2, numr, numi, l, &rmax);
00405     cmpadd_(sumr, sumi, numr, numi, sumr, sumi, l, &rmax);
00406     cnt += sigfig;
00407     ar += sigfig;
00408     cr += sigfig;
00409     goto L110;
00410 L190:
00411     arydiv_(sumr, sumi, denomr, denomi, &final, l, lnchf, &rmax, &bit);
00412      ret_val->r = final.r,  ret_val->i = final.i;
00413     return ;
00414 } /* chgf_ */
00415 
00416 /*     **************************************************************** */
00417 /*     *                                                              * */
00418 /*     *                 SUBROUTINE ARADD                             * */
00419 /*     *                                                              * */
00420 /*     *                                                              * */
00421 /*     *  Description : Accepts two arrays of numbers and returns     * */
00422 /*     *    the sum of the array.  Each array is holding the value    * */
00423 /*     *    of one number in the series.  The parameter L is the      * */
00424 /*     *    size of the array representing the number and RMAX is     * */
00425 /*     *    the actual number of digits needed to give the numbers    * */
00426 /*     *    the desired accuracy.                                     * */
00427 /*     *                                                              * */
00428 /*     *  Subprograms called: none                                    * */
00429 /*     *                                                              * */
00430 /*     **************************************************************** */
00431 /* Subroutine */ 
00432 int aradd_(const doublereal *a, const doublereal *b, 
00433         doublereal *c__, const integer *l, const doublereal *rmax)
00434 {
00435     /* System generated locals */
00436     integer i__1;
00437     doublereal d__1;
00438 
00439     /* Builtin functions */
00440 
00441     /* Local variables */
00442     static integer i__, j, ediff;
00443     static doublereal z__[779];
00444 
00445     /* Parameter adjustments */
00446     ++c__;
00447     ++b;
00448     ++a;
00449 
00450     /* Function Body */
00451     i__1 = *l + 1;
00452     for (i__ = 0; i__ <= i__1; ++i__) {
00453         z__[i__ + 1] = 0.;
00454 /* L110: */
00455     }
00456     d__1 = a[*l + 1] - b[*l + 1];
00457     ediff = (integer) d_nint(&d__1);
00458     if (abs(a[1]) < (float).5 || ediff <= -(*l)) {
00459         goto L111;
00460     }
00461     if (abs(b[1]) < (float).5 || ediff >= *l) {
00462         goto L113;
00463     }
00464     goto L115;
00465 L111:
00466     i__1 = *l + 1;
00467     for (i__ = -1; i__ <= i__1; ++i__) {
00468         c__[i__] = b[i__];
00469 /* L112: */
00470     }
00471     goto L311;
00472 L113:
00473     i__1 = *l + 1;
00474     for (i__ = -1; i__ <= i__1; ++i__) {
00475         c__[i__] = a[i__];
00476 /* L114: */
00477     }
00478     goto L311;
00479 L115:
00480     z__[0] = a[-1];
00481     if ((d__1 = a[-1] - b[-1], abs(d__1)) < (float).5) {
00482         goto L200;
00483     }
00484     if (ediff > 0) {
00485         z__[*l + 2] = a[*l + 1];
00486         goto L233;
00487     }
00488     if (ediff < 0) {
00489         z__[*l + 2] = b[*l + 1];
00490         z__[0] = b[-1];
00491         goto L266;
00492     }
00493     i__1 = *l;
00494     for (i__ = 1; i__ <= i__1; ++i__) {
00495         if (a[i__] > b[i__]) {
00496             z__[*l + 2] = a[*l + 1];
00497             goto L233;
00498         }
00499         if (a[i__] < b[i__]) {
00500             z__[*l + 2] = b[*l + 1];
00501             z__[0] = b[-1];
00502             goto L266;
00503         }
00504 /* L120: */
00505     }
00506     goto L300;
00507 L200:
00508     if (ediff > 0) {
00509         goto L203;
00510     }
00511     if (ediff < 0) {
00512         goto L207;
00513     }
00514     z__[*l + 2] = a[*l + 1];
00515     for (i__ = *l; i__ >= 1; --i__) {
00516         z__[i__ + 1] = a[i__] + b[i__] + z__[i__ + 1];
00517         if (z__[i__ + 1] >= *rmax) {
00518             z__[i__ + 1] -= *rmax;
00519             z__[i__] = 1.;
00520         }
00521 /* L201: */
00522     }
00523     if (z__[1] > (float).5) {
00524         for (i__ = *l; i__ >= 1; --i__) {
00525             z__[i__ + 1] = z__[i__];
00526 /* L202: */
00527         }
00528         z__[*l + 2] += 1.;
00529         z__[1] = 0.;
00530     }
00531     goto L300;
00532 L203:
00533     z__[*l + 2] = a[*l + 1];
00534     i__1 = ediff + 1;
00535     for (i__ = *l; i__ >= i__1; --i__) {
00536         z__[i__ + 1] = a[i__] + b[i__ - ediff] + z__[i__ + 1];
00537         if (z__[i__ + 1] >= *rmax) {
00538             z__[i__ + 1] -= *rmax;
00539             z__[i__] = 1.;
00540         }
00541 /* L204: */
00542     }
00543     for (i__ = ediff; i__ >= 1; --i__) {
00544         z__[i__ + 1] = a[i__] + z__[i__ + 1];
00545         if (z__[i__ + 1] >= *rmax) {
00546             z__[i__ + 1] -= *rmax;
00547             z__[i__] = 1.;
00548         }
00549 /* L205: */
00550     }
00551     if (z__[1] > (float).5) {
00552         for (i__ = *l; i__ >= 1; --i__) {
00553             z__[i__ + 1] = z__[i__];
00554 /* L206: */
00555         }
00556         z__[*l + 2] += 1;
00557         z__[1] = 0.;
00558     }
00559     goto L300;
00560 L207:
00561     z__[*l + 2] = b[*l + 1];
00562     i__1 = 1 - ediff;
00563     for (i__ = *l; i__ >= i__1; --i__) {
00564         z__[i__ + 1] = a[i__ + ediff] + b[i__] + z__[i__ + 1];
00565         if (z__[i__ + 1] >= *rmax) {
00566             z__[i__ + 1] -= *rmax;
00567             z__[i__] = 1.;
00568         }
00569 /* L208: */
00570     }
00571     for (i__ = -ediff; i__ >= 1; --i__) {
00572         z__[i__ + 1] = b[i__] + z__[i__ + 1];
00573         if (z__[i__ + 1] >= *rmax) {
00574             z__[i__ + 1] -= *rmax;
00575             z__[i__] = 1.;
00576         }
00577 /* L209: */
00578     }
00579     if (z__[1] > (float).5) {
00580         for (i__ = *l; i__ >= 1; --i__) {
00581             z__[i__ + 1] = z__[i__];
00582 /* L210: */
00583         }
00584         z__[*l + 2] += 1.;
00585         z__[1] = 0.;
00586     }
00587     goto L300;
00588 L233:
00589     if (ediff > 0) {
00590         goto L243;
00591     }
00592     for (i__ = *l; i__ >= 1; --i__) {
00593         z__[i__ + 1] = a[i__] - b[i__] + z__[i__ + 1];
00594         if (z__[i__ + 1] < 0.) {
00595             z__[i__ + 1] += *rmax;
00596             z__[i__] = -1.;
00597         }
00598 /* L234: */
00599     }
00600     goto L290;
00601 L243:
00602     i__1 = ediff + 1;
00603     for (i__ = *l; i__ >= i__1; --i__) {
00604         z__[i__ + 1] = a[i__] - b[i__ - ediff] + z__[i__ + 1];
00605         if (z__[i__ + 1] < 0.) {
00606             z__[i__ + 1] += *rmax;
00607             z__[i__] = -1.;
00608         }
00609 /* L244: */
00610     }
00611     for (i__ = ediff; i__ >= 1; --i__) {
00612         z__[i__ + 1] = a[i__] + z__[i__ + 1];
00613         if (z__[i__ + 1] < 0.) {
00614             z__[i__ + 1] += *rmax;
00615             z__[i__] = -1.;
00616         }
00617 /* L245: */
00618     }
00619     goto L290;
00620 L266:
00621     if (ediff < 0) {
00622         goto L276;
00623     }
00624     for (i__ = *l; i__ >= 1; --i__) {
00625         z__[i__ + 1] = b[i__] - a[i__] + z__[i__ + 1];
00626         if (z__[i__ + 1] < 0.) {
00627             z__[i__ + 1] += *rmax;
00628             z__[i__] = -1.;
00629         }
00630 /* L267: */
00631     }
00632     goto L290;
00633 L276:
00634     i__1 = 1 - ediff;
00635     for (i__ = *l; i__ >= i__1; --i__) {
00636         z__[i__ + 1] = b[i__] - a[i__ + ediff] + z__[i__ + 1];
00637         if (z__[i__ + 1] < 0.) {
00638             z__[i__ + 1] += *rmax;
00639             z__[i__] = -1.;
00640         }
00641 /* L277: */
00642     }
00643     for (i__ = -ediff; i__ >= 1; --i__) {
00644         z__[i__ + 1] = b[i__] + z__[i__ + 1];
00645         if (z__[i__ + 1] < 0.) {
00646             z__[i__ + 1] += *rmax;
00647             z__[i__] = -1.;
00648         }
00649 /* L278: */
00650     }
00651 L290:
00652     if (z__[2] > (float).5) {
00653         goto L300;
00654     }
00655     i__ = 1;
00656 L291:
00657     ++i__;
00658     if (z__[i__ + 1] < (float).5 && i__ < *l + 1) {
00659         goto L291;
00660     }
00661     if (i__ == *l + 1) {
00662         z__[0] = 1.;
00663         z__[*l + 2] = 0.;
00664         goto L300;
00665     }
00666 /* L292: */
00667     i__1 = *l + 1 - i__;
00668     for (j = 1; j <= i__1; ++j) {
00669         z__[j + 1] = z__[j + i__];
00670 /* L293: */
00671     }
00672     i__1 = *l;
00673     for (j = *l + 2 - i__; j <= i__1; ++j) {
00674         z__[j + 1] = 0.;
00675 /* L294: */
00676     }
00677     z__[*l + 2] = z__[*l + 2] - i__ + 1;
00678 L300:
00679     i__1 = *l + 1;
00680     for (i__ = -1; i__ <= i__1; ++i__) {
00681         c__[i__] = z__[i__ + 1];
00682 /* L310: */
00683     }
00684 L311:
00685     if (c__[1] < (float).5) {
00686         c__[-1] = 1.;
00687         c__[*l + 1] = 0.;
00688     }
00689     return 0;
00690 } /* aradd_ */
00691 
00692 /*     **************************************************************** */
00693 /*     *                                                              * */
00694 /*     *                 SUBROUTINE ARSUB                             * */
00695 /*     *                                                              * */
00696 /*     *                                                              * */
00697 /*     *  Description : Accepts two arrays and subtracts each element * */
00698 /*     *    in the second array from the element in the first array   * */
00699 /*     *    and returns the solution.  The parameters L and RMAX are  * */
00700 /*     *    the size of the array and the number of digits needed for * */
00701 /*     *    the accuracy, respectively.                               * */
00702 /*     *                                                              * */
00703 /*     *  Subprograms called: ARADD                                   * */
00704 /*     *                                                              * */
00705 /*     **************************************************************** */
00706 /* Subroutine */ 
00707 int arsub_(const doublereal *a, const doublereal *b, 
00708         doublereal *c__, const integer *l, const doublereal *rmax)
00709 {
00710     /* System generated locals */
00711     integer i__1;
00712 
00713     /* Local variables */
00714     static integer i__;
00715     static doublereal b2[779];
00716 
00717     /* Parameter adjustments */
00718     ++c__;
00719     ++b;
00720     ++a;
00721 
00722     /* Function Body */
00723     i__1 = *l + 1;
00724     for (i__ = -1; i__ <= i__1; ++i__) {
00725         b2[i__ + 1] = b[i__];
00726 /* L100: */
00727     }
00728     b2[0] *= -1.;
00729     aradd_(&a[-1], b2, &c__[-1], l, rmax);
00730     return 0;
00731 } /* arsub_ */
00732 
00733 /*     **************************************************************** */
00734 /*     *                                                              * */
00735 /*     *                 SUBROUTINE ARMULT                            * */
00736 /*     *                                                              * */
00737 /*     *                                                              * */
00738 /*     *  Description : Accepts two arrays and returns the product.   * */
00739 /*     *    L and RMAX are the size of the arrays and the number of   * */
00740 /*     *    digits needed to represent the numbers with the required  * */
00741 /*     *    accuracy.                                                 * */
00742 /*     *                                                              * */
00743 /*     *  Subprograms called: none                                    * */
00744 /*     *                                                              * */
00745 /*     **************************************************************** */
00746 /* Subroutine */ int armult_(const doublereal *a, const doublereal *b, doublereal *c__, const integer *l, const doublereal *rmax)
00747 {
00748     /* System generated locals */
00749     integer i__1;
00750     doublereal d__1;
00751 
00752     /* Builtin functions */
00753 
00754     /* Local variables */
00755     static doublereal rmax2;
00756     static integer i__;
00757     static doublereal z__[779], carry, b2;
00758 
00759     /* Parameter adjustments */
00760     ++c__;
00761     ++a;
00762 
00763     /* Function Body */
00764     rmax2 = 1. / *rmax;
00765     z__[0] = d_sign(&c_b53, b) * a[-1];
00766     b2 = abs(*b);
00767     z__[*l + 2] = a[*l + 1];
00768     i__1 = *l;
00769     for (i__ = 0; i__ <= i__1; ++i__) {
00770         z__[i__ + 1] = 0.;
00771 /* L100: */
00772     }
00773     if (b2 <= 1e-10 || a[1] <= 1e-10) {
00774         z__[0] = 1.;
00775         z__[*l + 2] = 0.;
00776         goto L198;
00777     }
00778     for (i__ = *l; i__ >= 1; --i__) {
00779         z__[i__ + 1] = a[i__] * b2 + z__[i__ + 1];
00780         if (z__[i__ + 1] >= *rmax) {
00781             d__1 = z__[i__ + 1] / *rmax;
00782             carry = d_int(&d__1);
00783             z__[i__ + 1] -= carry * *rmax;
00784             z__[i__] = carry;
00785         }
00786 /* L110: */
00787     }
00788     if (z__[1] < (float).5) {
00789         goto L150;
00790     }
00791     for (i__ = *l; i__ >= 1; --i__) {
00792         z__[i__ + 1] = z__[i__];
00793 /* L120: */
00794     }
00795     z__[*l + 2] += 1.;
00796     z__[1] = 0.;
00797 L150:
00798 L198:
00799     i__1 = *l + 1;
00800     for (i__ = -1; i__ <= i__1; ++i__) {
00801         c__[i__] = z__[i__ + 1];
00802 /* L199: */
00803     }
00804     if (c__[1] < (float).5) {
00805         c__[-1] = 1.;
00806         c__[*l + 1] = 0.;
00807     }
00808     return 0;
00809 } /* armult_ */
00810 
00811 /*     **************************************************************** */
00812 /*     *                                                              * */
00813 /*     *                 SUBROUTINE CMPADD                            * */
00814 /*     *                                                              * */
00815 /*     *                                                              * */
00816 /*     *  Description : Takes two arrays representing one real and    * */
00817 /*     *    one imaginary part, and adds two arrays representing      * */
00818 /*     *    another complex number and returns two array holding the  * */
00819 /*     *    complex sum.                                              * */
00820 /*     *              (CR,CI) = (AR+BR, AI+BI)                        * */
00821 /*     *                                                              * */
00822 /*     *  Subprograms called: ARADD                                   * */
00823 /*     *                                                              * */
00824 /*     **************************************************************** */
00825 /* 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)
00826 {
00827 
00828     /* Parameter adjustments */
00829     ++ci;
00830     ++cr;
00831     ++bi;
00832     ++br;
00833     ++ai;
00834     ++ar;
00835 
00836     /* Function Body */
00837     aradd_(&ar[-1], &br[-1], &cr[-1], l, rmax);
00838     aradd_(&ai[-1], &bi[-1], &ci[-1], l, rmax);
00839     return 0;
00840 } /* cmpadd_ */
00841 
00842 /*     **************************************************************** */
00843 /*     *                                                              * */
00844 /*     *                 SUBROUTINE CMPSUB                            * */
00845 /*     *                                                              * */
00846 /*     *                                                              * */
00847 /*     *  Description : Takes two arrays representing one real and    * */
00848 /*     *    one imaginary part, and subtracts two arrays representing * */
00849 /*     *    another complex number and returns two array holding the  * */
00850 /*     *    complex sum.                                              * */
00851 /*     *              (CR,CI) = (AR+BR, AI+BI)                        * */
00852 /*     *                                                              * */
00853 /*     *  Subprograms called: ARADD                                   * */
00854 /*     *                                                              * */
00855 /*     **************************************************************** */
00856 /* Subroutine */ int cmpsub_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci, integer *l, doublereal *rmax)
00857 {
00858     /* Parameter adjustments */
00859     ++ci;
00860     ++cr;
00861     ++bi;
00862     ++br;
00863     ++ai;
00864     ++ar;
00865 
00866     /* Function Body */
00867     arsub_(&ar[-1], &br[-1], &cr[-1], l, rmax);
00868     arsub_(&ai[-1], &bi[-1], &ci[-1], l, rmax);
00869     return 0;
00870 } /* cmpsub_ */
00871 
00872 /*     **************************************************************** */
00873 /*     *                                                              * */
00874 /*     *                 SUBROUTINE CMPMUL                            * */
00875 /*     *                                                              * */
00876 /*     *                                                              * */
00877 /*     *  Description : Takes two arrays representing one real and    * */
00878 /*     *    one imaginary part, and multiplies it with two arrays     * */
00879 /*     *    representing another complex number and returns the       * */
00880 /*     *    complex product.                                          * */
00881 /*     *                                                              * */
00882 /*     *  Subprograms called: ARMULT, ARSUB, ARADD                    * */
00883 /*     *                                                              * */
00884 /*     **************************************************************** */
00885 /* 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)
00886 {
00887     static doublereal d1[779], d2[779];
00888 
00889     /* Parameter adjustments */
00890     ++ci;
00891     ++cr;
00892     ++ai;
00893     ++ar;
00894 
00895     /* Function Body */
00896     armult_(&ar[-1], br, d1, l, rmax);
00897     armult_(&ai[-1], bi, d2, l, rmax);
00898     arsub_(d1, d2, &cr[-1], l, rmax);
00899     armult_(&ar[-1], bi, d1, l, rmax);
00900     armult_(&ai[-1], br, d2, l, rmax);
00901     aradd_(d1, d2, &ci[-1], l, rmax);
00902     return 0;
00903 } /* cmpmul_ */
00904 
00905 /*     **************************************************************** */
00906 /*     *                                                              * */
00907 /*     *                 SUBROUTINE ARYDIV                            * */
00908 /*     *                                                              * */
00909 /*     *                                                              * */
00910 /*     *  Description : Returns the double precision complex number   * */
00911 /*     *    resulting from the division of four arrays, representing  * */
00912 /*     *    two complex numbers.  The number returned will be in one  * */
00913 /*     *    two different forms.  Either standard scientific or as    * */
00914 /*     *    the log of the number.                                    * */
00915 /*     *                                                              * */
00916 /*     *  Subprograms called: CONV21, CONV12, EADD, ECPDIV, EMULT     * */
00917 /*     *                                                              * */
00918 /*     **************************************************************** */
00919 /* 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)
00920 {
00921     /* System generated locals */
00922     doublereal d__1;
00923     doublecomplex z__1;
00924 
00925     /* Builtin functions */
00926 
00927     /* Local variables */
00928     static integer rexp;
00929     static doublereal x;
00930     static doublereal e1, e2, e3;
00931     static doublereal n1, n2, n3, x1, x2, ae[4] /* was [2][2] */, be[4] /* 
00932             was [2][2] */, ce[4]        /* was [2][2] */;
00933     static integer ii10, ir10;
00934     static doublereal ri10, phi, rr10, dum1, dum2;
00935     /* Parameter adjustments */
00936     ++bi;
00937     ++br;
00938     ++ai;
00939     ++ar;
00940 
00941     /* Function Body */
00942     rexp = *bit / 2;
00943     x = rexp * (ar[*l + 1] - 2);
00944     rr10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
00945     ir10 = (integer) rr10;
00946     rr10 -= ir10;
00947     x = rexp * (ai[*l + 1] - 2);
00948     ri10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
00949     ii10 = (integer) ri10;
00950     ri10 -= ii10;
00951     d__1 = ar[1] * *rmax * *rmax + ar[2] * *rmax + ar[3];
00952     dum1 = d_sign(&d__1, &ar[-1]);
00953     d__1 = ai[1] * *rmax * *rmax + ai[2] * *rmax + ai[3];
00954     dum2 = d_sign(&d__1, &ai[-1]);
00955     dum1 *= pow_dd(&c_b65, &rr10);
00956     dum2 *= pow_dd(&c_b65, &ri10);
00957     z__1.r = dum1, z__1.i = dum2;
00958     conv12_(&z__1, ae);
00959     ae[2] += ir10;
00960     ae[3] += ii10;
00961     x = rexp * (br[*l + 1] - 2);
00962     rr10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
00963     ir10 = (integer) rr10;
00964     rr10 -= ir10;
00965     x = rexp * (bi[*l + 1] - 2);
00966     ri10 = x * d_lg10(&c_b8) / d_lg10(&c_b65);
00967     ii10 = (integer) ri10;
00968     ri10 -= ii10;
00969     d__1 = br[1] * *rmax * *rmax + br[2] * *rmax + br[3];
00970     dum1 = d_sign(&d__1, &br[-1]);
00971     d__1 = bi[1] * *rmax * *rmax + bi[2] * *rmax + bi[3];
00972     dum2 = d_sign(&d__1, &bi[-1]);
00973     dum1 *= pow_dd(&c_b65, &rr10);
00974     dum2 *= pow_dd(&c_b65, &ri10);
00975     z__1.r = dum1, z__1.i = dum2;
00976     conv12_(&z__1, be);
00977     be[2] += ir10;
00978     be[3] += ii10;
00979     ecpdiv_(ae, be, ce);
00980     if (*lnchf == 0) {
00981         conv21_(ce, c__);
00982     } else {
00983         emult_(ce, &ce[2], ce, &ce[2], &n1, &e1);
00984         emult_(&ce[1], &ce[3], &ce[1], &ce[3], &n2, &e2);
00985         eadd_(&n1, &e1, &n2, &e2, &n3, &e3);
00986         n1 = ce[0];
00987         e1 = ce[2] - ce[3];
00988         x2 = ce[1];
00989         if (e1 > 74.) {
00990             x1 = 1e75;
00991         } else if (e1 < -74.) {
00992             x1 = 0.;
00993         } else {
00994             x1 = n1 * pow_dd(&c_b65, &e1);
00995         }
00996         phi = atan2(x2, x1);
00997         d__1 = (log(n3) + e3 * log(10.)) * .5;
00998         z__1.r = d__1, z__1.i = phi;
00999         c__->r = z__1.r, c__->i = z__1.i;
01000     }
01001     return 0;
01002 } /* arydiv_ */
01003 
01004 /*     **************************************************************** */
01005 /*     *                                                              * */
01006 /*     *                 SUBROUTINE EMULT                             * */
01007 /*     *                                                              * */
01008 /*     *                                                              * */
01009 /*     *  Description : Takes one base and exponent and multiplies it * */
01010 /*     *    by another numbers base and exponent to give the product  * */
01011 /*     *    in the form of base and exponent.                         * */
01012 /*     *                                                              * */
01013 /*     *  Subprograms called: none                                    * */
01014 /*     *                                                              * */
01015 /*     **************************************************************** */
01016 /* Subroutine */ int emult_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
01017 {
01018     *nf = *n1 * *n2;
01019     *ef = *e1 + *e2;
01020     if (abs(*nf) >= 10.) {
01021         *nf /= 10.;
01022         *ef += 1.;
01023     }
01024     return 0;
01025 } /* emult_ */
01026 
01027 /*     **************************************************************** */
01028 /*     *                                                              * */
01029 /*     *                 SUBROUTINE EDIV                              * */
01030 /*     *                                                              * */
01031 /*     *                                                              * */
01032 /*     *  Description : returns the solution in the form of base and  * */
01033 /*     *    exponent of the division of two exponential numbers.      * */
01034 /*     *                                                              * */
01035 /*     *  Subprograms called: none                                    * */
01036 /*     *                                                              * */
01037 /*     **************************************************************** */
01038 /* Subroutine */ int ediv_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
01039 {
01040     *nf = *n1 / *n2;
01041     *ef = *e1 - *e2;
01042     if (abs(*nf) < 1. && *nf != 0.) {
01043         *nf *= 10.;
01044         *ef += -1.;
01045     }
01046     return 0;
01047 } /* ediv_ */
01048 
01049 /*     **************************************************************** */
01050 /*     *                                                              * */
01051 /*     *                 SUBROUTINE EADD                              * */
01052 /*     *                                                              * */
01053 /*     *                                                              * */
01054 /*     *  Description : Returns the sum of two numbers in the form    * */
01055 /*     *    of a base and an exponent.                                * */
01056 /*     *                                                              * */
01057 /*     *  Subprograms called: none                                    * */
01058 /*     *                                                              * */
01059 /*     **************************************************************** */
01060 /* Subroutine */ int eadd_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
01061 {
01062     /* Builtin functions */
01063 
01064     /* Local variables */
01065     static doublereal ediff;
01066 
01067     ediff = *e1 - *e2;
01068     if (ediff > 36.) {
01069         *nf = *n1;
01070         *ef = *e1;
01071     } else if (ediff < -36.) {
01072         *nf = *n2;
01073         *ef = *e2;
01074     } else {
01075         *nf = *n1 * pow_dd(&c_b65, &ediff) + *n2;
01076         *ef = *e2;
01077 L400:
01078         if (abs(*nf) < 10.) {
01079             goto L410;
01080         }
01081         *nf /= 10.;
01082         *ef += 1.;
01083         goto L400;
01084 L410:
01085         if (abs(*nf) >= 1. || *nf == 0.) {
01086             goto L420;
01087         }
01088         *nf *= 10.;
01089         *ef += -1.;
01090         goto L410;
01091     }
01092 L420:
01093     return 0;
01094 } /* eadd_ */
01095 
01096 /*     **************************************************************** */
01097 /*     *                                                              * */
01098 /*     *                 SUBROUTINE ESUB                              * */
01099 /*     *                                                              * */
01100 /*     *                                                              * */
01101 /*     *  Description : Returns the solution to the subtraction of    * */
01102 /*     *    two numbers in the form of base and exponent.             * */
01103 /*     *                                                              * */
01104 /*     *  Subprograms called: EADD                                    * */
01105 /*     *                                                              * */
01106 /*     **************************************************************** */
01107 /* Subroutine */ int esub_(doublereal *n1, doublereal *e1, doublereal *n2, doublereal *e2, doublereal *nf, doublereal *ef)
01108 {
01109     /* System generated locals */
01110     doublereal d__1;
01111 
01112     /* Local variables */
01113 
01114     d__1 = *n2 * -1.;
01115     eadd_(n1, e1, &d__1, e2, nf, ef);
01116     return 0;
01117 } /* esub_ */
01118 
01119 /*     **************************************************************** */
01120 /*     *                                                              * */
01121 /*     *                 SUBROUTINE CONV12                            * */
01122 /*     *                                                              * */
01123 /*     *                                                              * */
01124 /*     *  Description : Converts a number from complex notation to a  * */
01125 /*     *    form of a 2x2 real array.                                 * */
01126 /*     *                                                              * */
01127 /*     *  Subprograms called: none                                    * */
01128 /*     *                                                              * */
01129 /*     **************************************************************** */
01130 /* Subroutine */ int conv12_(doublecomplex *cn, doublereal *cae)
01131 {
01132     /* Builtin functions */
01133 
01134     /* Parameter adjustments */
01135     cae -= 3;
01136 
01137     /* Function Body */
01138     cae[3] = cn->r;
01139     cae[5] = 0.;
01140 L300:
01141     if (abs(cae[3]) < 10.) {
01142         goto L310;
01143     }
01144     cae[3] /= 10.;
01145     cae[5] += 1.;
01146     goto L300;
01147 L310:
01148     if (abs(cae[3]) >= 1. || cae[3] == 0.) {
01149         goto L320;
01150     }
01151     cae[3] *= 10.;
01152     cae[5] += -1.;
01153     goto L310;
01154 L320:
01155     cae[4] = d_imag(cn);
01156     cae[6] = 0.;
01157 L330:
01158     if (abs(cae[4]) < 10.) {
01159         goto L340;
01160     }
01161     cae[4] /= 10.;
01162     cae[6] += 1.;
01163     goto L330;
01164 L340:
01165     if (abs(cae[4]) >= 1. || cae[4] == 0.) {
01166         goto L350;
01167     }
01168     cae[4] *= 10.;
01169     cae[6] += -1.;
01170     goto L340;
01171 L350:
01172     return 0;
01173 } /* conv12_ */
01174 
01175 /*     **************************************************************** */
01176 /*     *                                                              * */
01177 /*     *                 SUBROUTINE CONV21                            * */
01178 /*     *                                                              * */
01179 /*     *                                                              * */
01180 /*     *  Description : Converts a number represented in a 2x2 real   * */
01181 /*     *    array to the form of a complex number.                    * */
01182 /*     *                                                              * */
01183 /*     *  Subprograms called: none                                    * */
01184 /*     *                                                              * */
01185 /*     **************************************************************** */
01186 /* Subroutine */ int conv21_(doublereal *cae, doublecomplex *cn)
01187 {
01188     /* System generated locals */
01189     doublereal d__1, d__2;
01190     doublecomplex z__1;
01191 
01192     /* Builtin functions */
01193 
01194     /* Parameter adjustments */
01195     cae -= 3;
01196 
01197     /* Function Body */
01198     if (cae[5] > 75. || cae[6] > 75.) {
01199         cn->r = 1e75, cn->i = 1e75;
01200     } else if (cae[6] < -75.) {
01201         d__1 = cae[3] * pow_dd(&c_b65, &cae[5]);
01202         z__1.r = d__1, z__1.i = 0.;
01203         cn->r = z__1.r, cn->i = z__1.i;
01204     } else {
01205         d__1 = cae[3] * pow_dd(&c_b65, &cae[5]);
01206         d__2 = cae[4] * pow_dd(&c_b65, &cae[6]);
01207         z__1.r = d__1, z__1.i = d__2;
01208         cn->r = z__1.r, cn->i = z__1.i;
01209     }
01210     return 0;
01211 } /* conv21_ */
01212 
01213 /*     **************************************************************** */
01214 /*     *                                                              * */
01215 /*     *                 SUBROUTINE ECPMUL                            * */
01216 /*     *                                                              * */
01217 /*     *                                                              * */
01218 /*     *  Description : Multiplies two numbers which are each         * */
01219 /*     *    represented in the form of a two by two array and returns * */
01220 /*     *    the solution in the same form.                            * */
01221 /*     *                                                              * */
01222 /*     *  Subprograms called: EMULT, ESUB, EADD                       * */
01223 /*     *                                                              * */
01224 /*     **************************************************************** */
01225 /* Subroutine */ int ecpmul_(doublereal *a, doublereal *b, doublereal *c__)
01226 {
01227     static doublereal c2[4]     /* was [2][2] */, e1, e2;
01228     static doublereal n1, n2;
01229     /* Parameter adjustments */
01230     c__ -= 3;
01231     b -= 3;
01232     a -= 3;
01233 
01234     /* Function Body */
01235     emult_(&a[3], &a[5], &b[3], &b[5], &n1, &e1);
01236     emult_(&a[4], &a[6], &b[4], &b[6], &n2, &e2);
01237     esub_(&n1, &e1, &n2, &e2, c2, &c2[2]);
01238     emult_(&a[3], &a[5], &b[4], &b[6], &n1, &e1);
01239     emult_(&a[4], &a[6], &b[3], &b[5], &n2, &e2);
01240     eadd_(&n1, &e1, &n2, &e2, &c__[4], &c__[6]);
01241     c__[3] = c2[0];
01242     c__[5] = c2[2];
01243     return 0;
01244 } /* ecpmul_ */
01245 
01246 /*     **************************************************************** */
01247 /*     *                                                              * */
01248 /*     *                 SUBROUTINE ECPDIV                            * */
01249 /*     *                                                              * */
01250 /*     *                                                              * */
01251 /*     *  Description : Divides two numbers and returns the solution. * */
01252 /*     *    All numbers are represented by a 2x2 array.               * */
01253 /*     *                                                              * */
01254 /*     *  Subprograms called: EADD, ECPMUL, EDIV, EMULT               * */
01255 /*     *                                                              * */
01256 /*     **************************************************************** */
01257 /* Subroutine */ int ecpdiv_(doublereal *a, doublereal *b, doublereal *c__)
01258 {
01259     static doublereal b2[4]     /* was [2][2] */, c2[4] /* was [2][2] */, e1, 
01260             e2, e3;
01261     static doublereal n1, n2, n3;
01262     /* Parameter adjustments */
01263     c__ -= 3;
01264     b -= 3;
01265     a -= 3;
01266 
01267     /* Function Body */
01268     b2[0] = b[3];
01269     b2[2] = b[5];
01270     b2[1] = b[4] * -1.;
01271     b2[3] = b[6];
01272     ecpmul_(&a[3], b2, c2);
01273     emult_(&b[3], &b[5], &b[3], &b[5], &n1, &e1);
01274     emult_(&b[4], &b[6], &b[4], &b[6], &n2, &e2);
01275     eadd_(&n1, &e1, &n2, &e2, &n3, &e3);
01276     ediv_(c2, &c2[2], &n3, &e3, &c__[3], &c__[5]);
01277     ediv_(&c2[1], &c2[3], &n3, &e3, &c__[4], &c__[6]);
01278     return 0;
01279 } /* ecpdiv_ */
01280 
01281 /* Main program */ int MAIN__(void)
01282 {
01283     /* Format strings */
01284     static char fmt_300[] = "(1x,\002 RESULT FROM CONHYP=\002,1p,2d25.12)";
01285     static char fmt_301[] = "(1x,\002 EXPECTED RESULT=\002,3x,1p,2d25.12)";
01286 
01287     /* System generated locals */
01288     doublecomplex z__1;
01289     olist o__1;
01290 
01291     /* Builtin functions */
01292 
01293     /* Local variables */
01294     static doublecomplex a, b, z__;
01295     static integer lnchf;
01296     static doublereal ai, ar;
01297     static integer ip;
01298     static doublecomplex chf;
01299 
01300     /* Fortran I/O blocks */
01301     static cilist io___88 = { 0, 5, 0, 0, 0 };
01302     static cilist io___90 = { 0, 6, 0, 0, 0 };
01303     static cilist io___91 = { 0, 5, 0, 0, 0 };
01304     static cilist io___93 = { 0, 6, 0, 0, 0 };
01305     static cilist io___94 = { 0, 5, 0, 0, 0 };
01306     static cilist io___96 = { 0, 6, 0, 0, 0 };
01307     static cilist io___97 = { 0, 5, 0, 0, 0 };
01308     static cilist io___99 = { 0, 6, 0, 0, 0 };
01309     static cilist io___100 = { 0, 5, 0, 0, 0 };
01310     static cilist io___102 = { 0, 6, 0, 0, 0 };
01311     static cilist io___104 = { 0, 6, 0, fmt_300, 0 };
01312     static cilist io___107 = { 0, 6, 0, fmt_301, 0 };
01313 
01314 
01315 
01316     o__1.oerr = 0;
01317     o__1.ounit = 5;
01318     o__1.ofnmlen = 4;
01319     o__1.ofnm = "DATA";
01320     o__1.orl = 0;
01321     o__1.osta = "OLD";
01322     o__1.oacc = 0;
01323     o__1.ofm = 0;
01324     o__1.oblnk = 0;
01325     f_open(&o__1);
01326     s_rsle(&io___88);
01327     do_lio(&c__7, &c__1, (char *)&a, (ftnlen)sizeof(doublecomplex));
01328     e_rsle();
01329     s_wsle(&io___90);
01330     do_lio(&c__9, &c__1, " A= ", (ftnlen)4);
01331     do_lio(&c__7, &c__1, (char *)&a, (ftnlen)sizeof(doublecomplex));
01332     e_wsle();
01333     s_rsle(&io___91);
01334     do_lio(&c__7, &c__1, (char *)&b, (ftnlen)sizeof(doublecomplex));
01335     e_rsle();
01336     s_wsle(&io___93);
01337     do_lio(&c__9, &c__1, " B= ", (ftnlen)4);
01338     do_lio(&c__7, &c__1, (char *)&b, (ftnlen)sizeof(doublecomplex));
01339     e_wsle();
01340     s_rsle(&io___94);
01341     do_lio(&c__7, &c__1, (char *)&z__, (ftnlen)sizeof(doublecomplex));
01342     e_rsle();
01343     s_wsle(&io___96);
01344     do_lio(&c__9, &c__1, " Z= ", (ftnlen)4);
01345     do_lio(&c__7, &c__1, (char *)&z__, (ftnlen)sizeof(doublecomplex));
01346     e_wsle();
01347     s_rsle(&io___97);
01348     do_lio(&c__3, &c__1, (char *)&lnchf, (ftnlen)sizeof(integer));
01349     e_rsle();
01350     s_wsle(&io___99);
01351     do_lio(&c__9, &c__1, " LNCHF= ", (ftnlen)8);
01352     do_lio(&c__3, &c__1, (char *)&lnchf, (ftnlen)sizeof(integer));
01353     e_wsle();
01354     s_rsle(&io___100);
01355     do_lio(&c__3, &c__1, (char *)&ip, (ftnlen)sizeof(integer));
01356     e_rsle();
01357     s_wsle(&io___102);
01358     do_lio(&c__9, &c__1, " IP= ", (ftnlen)5);
01359     do_lio(&c__3, &c__1, (char *)&ip, (ftnlen)sizeof(integer));
01360     e_wsle();
01361     conhyp_(&z__1, &a, &b, &z__, &lnchf, &ip);
01362     chf.r = z__1.r, chf.i = z__1.i;
01363     s_wsfe(&io___104);
01364     do_fio(&c__2, (char *)&chf, (ftnlen)sizeof(doublereal));
01365     e_wsfe();
01366     ar = 2.31145634403113e-12;
01367     ai = -1.96169649634905e-11;
01368     s_wsfe(&io___107);
01369     do_fio(&c__1, (char *)&ar, (ftnlen)sizeof(doublereal));
01370     do_fio(&c__1, (char *)&ai, (ftnlen)sizeof(doublereal));
01371     e_wsfe();
01372     s_stop("", (ftnlen)0);
01373     return 0;
01374 } /* MAIN__ */
01375 
01376 /* Main program alias */ int sample_ (void) { return MAIN__ (); }
01377 #ifdef __cplusplus
01378 }
01379 #endif
01380