TOMS404.C

Go to the documentation of this file.

/* toms404.f -- translated by f2c (version 19970219).
   You must link the resulting object file with the libraries:
        -lf2c -lm   (in that order)
*/

/* $Id: TOMS404.C,v 1.2.2.12 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"
/* Table of constant values */

static integer c__2 = 2;

/*     ALGORITHM 404 COLLECTED ALGORITHMS FROM ACM. */
/*     ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 01, */
/*     P. 048. */
/* Complex */ /*VOID*/
void cgamma_(complex *ret_val, complex *z__)
{
    /* Format strings */
    static char fmt_900[] = "(1x,2e14.7,10x,\002ARGUMENT OF GAMMA FUNCTION I\
S TOO CLOSE TO A POLE\002)";
    static char fmt_910[] = "(\002 ERROR - STIRLING*S SERIES HAS NOT CONVERG\
ED\002/14x,\002Z = \002,2e14.7)";

    /* System generated locals */
    integer i__1;
    real r__1;
    complex q__1, q__2, q__3, q__4, q__5, q__6;

    /* Builtin functions: declared in prototypes.h */

    /* Local variables */
    static complex term;
    static integer iout;
    static complex a;
    static real c__[12];
    static integer i__, j, m;
    static complex t;
    static real x, y, xdist;
    static complex pi, zm, tt;
    static logical reflek;
    static complex den;
    static real tol;
    static complex sum;

    /* Fortran I/O blocks */
    static cilist io___9 = { 0, 0, 0, fmt_900, 0 };
    static cilist io___18 = { 0, 0, 0, fmt_910, 0 };


/* SET IOUT FOR PROPER OUTPUT CHANNEL OF COMPUTER SYSTEM FOR */
/* ERROR MESSAGES */
    iout = 3;
    pi.r = (float)3.14159265359, pi.i = (float)0.;
    x = z__->r;
    y = r_imag(z__);
/* TOL = LIMIT OF PRECISION OF COMPUTER SYSTEM IN SINGLE PRECISI */
    tol = (float)1e-7;
    reflek = TRUE_;
/* DETERMINE WHETHER Z IS TOO CLOSE TO A POLE */
/* CHECK WHETHER TOO CLOSE TO ORIGIN */
    if (x >= tol) {
        goto L20;
    }
/* FIND THE NEAREST POLE AND COMPUTE DISTANCE TO IT */
    xdist = x - (integer) (x - (float).5);
    q__1.r = xdist, q__1.i = y;
    zm.r = q__1.r, zm.i = q__1.i;
    if (c_abs(&zm) >= tol) {
        goto L10;
    }
/* IF Z IS TOO CLOSE TO A POLE, PRINT ERROR MESSAGE AND RETURN */
/* WITH CGAMMA = (1.E7,0.0E0) */
    io___9.ciunit = iout;
    s_wsfe(&io___9);
    do_fio(&c__2, (char *)&(*z__), (ftnlen)sizeof(real));
    e_wsfe();
     ret_val->r = (float)1e7,  ret_val->i = (float)0.;
    return ;
/* FOR REAL(Z) NEGATIVE EMPLOY THE REFLECTION FORMULA */
/* GAMMA(Z) = PI/(SIN(PI*Z)*GAMMA(1-Z)) */
/* AND COMPUTE GAMMA(1-Z).  NOTE REFLEK IS A TAG TO INDICATE THA */
/* THIS RELATION MUST BE USED LATER. */
L10:
    if (x >= (float)0.) {
        goto L20;
    }
    reflek = FALSE_;
    q__1.r = (float)1. - z__->r, q__1.i = (float)0. - z__->i;
    z__->r = q__1.r, z__->i = q__1.i;
    x = (float)1. - x;
    y = -y;
/* IF Z IS NOT TOO CLOSE TO A POLE, MAKE REAL(Z)>10 AND ARG(Z)<P */
L20:
    m = 0;
L40:
    if (x >= (float)10.) {
        goto L50;
    }
    x += (float)1.;
    ++m;
    goto L40;
L50:
    if (dabs(y) < x) {
        goto L60;
    }
    x += (float)1.;
    ++m;
    goto L50;
L60:
    q__1.r = x, q__1.i = y;
    t.r = q__1.r, t.i = q__1.i;
    q__1.r = t.r * t.r - t.i * t.i, q__1.i = t.r * t.i + t.i * t.r;
    tt.r = q__1.r, tt.i = q__1.i;
    den.r = t.r, den.i = t.i;
/* COEFFICIENTS IN STIRLING*S APPROXIMATION FOR LN(GAMMA(T)) */
    c__[0] = (float).08333333333333333;
    c__[1] = (float)-.002777777777777778;
    c__[2] = (float)7.936507936507937e-4;
    c__[3] = (float)-5.952380952380953e-4;
    c__[4] = (float)8.417508417508417e-4;
    c__[5] = (float)-.001917526917526918;
    c__[6] = (float).00641025641025641;
    c__[7] = (float)-.02955065359477124;
    c__[8] = (float).1796443723688306;
    c__[9] = (float)-1.392432216905901;
    c__[10] = (float)13.40286404416839;
    q__4.r = t.r - (float).5, q__4.i = t.i + (float)0.;
    c_log(&q__5, &t);
    q__3.r = q__4.r * q__5.r - q__4.i * q__5.i, q__3.i = q__4.r * q__5.i +
            q__4.i * q__5.r;
    q__2.r = q__3.r - t.r, q__2.i = q__3.i - t.i;
    r__1 = log((float)6.28318) * (float).5;
    q__6.r = r__1, q__6.i = (float)0.;
    q__1.r = q__2.r + q__6.r, q__1.i = q__2.i + q__6.i;
    sum.r = q__1.r, sum.i = q__1.i;
    j = 1;
L70:
    i__1 = j - 1;
    q__2.r = c__[i__1], q__2.i = (float)0.;
    c_div(&q__1, &q__2, &den);
    term.r = q__1.r, term.i = q__1.i;
/* TEST REAL AND IMAGINARY PARTS OF LN(GAMMA(Z)) SEPARATELY FOR */
/* CONVERGENCE.  IF Z IS REAL SKIP IMAGINARY PART OF CHECK. */
    if ((r__1 = term.r / sum.r, dabs(r__1)) >= tol) {
        goto L80;
    }
    if (y == (float)0.) {
        goto L100;
    }
    if ((r__1 = r_imag(&term) / r_imag(&sum), dabs(r__1)) < tol) {
        goto L100;
    }
L80:
    q__1.r = sum.r + term.r, q__1.i = sum.i + term.i;
    sum.r = q__1.r, sum.i = q__1.i;
    ++j;
    q__1.r = den.r * tt.r - den.i * tt.i, q__1.i = den.r * tt.i + den.i *
            tt.r;
    den.r = q__1.r, den.i = q__1.i;
/* TEST FOR NONCONVERGENCE */
    if (j == 12) {
        goto L90;
    }
    goto L70;
/* STIRLING*S SERIES DID NOT CONVERGE.  PRINT ERROR MESSAGE AND */
/* PROCEDE. */
L90:
    io___18.ciunit = iout;
    s_wsfe(&io___18);
    do_fio(&c__2, (char *)&(*z__), (ftnlen)sizeof(real));
    e_wsfe();
/* RECURSION RELATION USED TO OBTAIN LN(GAMMA(Z)) */
/* LN(GAMMA(Z)) = LN(GAMMA(Z+M)/(Z*(Z+1)*...*(Z+M-1))) */
/* = LN(GAMMA(Z+M)-LN(Z)-LN(Z+1)-...-LN(Z+M */
L100:
    if (m == 0) {
        goto L120;
    }
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        r__1 = i__ * (float)1. - (float)1.;
        q__1.r = r__1, q__1.i = (float)0.;
        a.r = q__1.r, a.i = q__1.i;
/* L110: */
        q__3.r = z__->r + a.r, q__3.i = z__->i + a.i;
        c_log(&q__2, &q__3);
        q__1.r = sum.r - q__2.r, q__1.i = sum.i - q__2.i;
        sum.r = q__1.r, sum.i = q__1.i;
    }
/* CHECK TO SEE IF REFLECTION FORMULA SHOULD BE USED */
L120:
    if (reflek) {
        goto L130;
    }
    q__5.r = pi.r * z__->r - pi.i * z__->i, q__5.i = pi.r * z__->i + pi.i *
            z__->r;
    c_sin(&q__4, &q__5);
    c_div(&q__3, &pi, &q__4);
    c_log(&q__2, &q__3);
    q__1.r = q__2.r - sum.r, q__1.i = q__2.i - sum.i;
    sum.r = q__1.r, sum.i = q__1.i;
    q__1.r = (float)1. - z__->r, q__1.i = (float)0. - z__->i;
    z__->r = q__1.r, z__->i = q__1.i;
L130:
    c_exp(&q__1, &sum);
     ret_val->r = q__1.r,  ret_val->i = q__1.i;
    return ;
} /* cgamma_ */

#ifdef __cplusplus
}
#endif