TBCI Numerical high perf. C++ Library: TOMS404.C Source File

00001 
00004 /* toms404.f -- translated by f2c (version 19970219).
00005    You must link the resulting object file with the libraries:
00006         -lf2c -lm   (in that order)
00007 */
00008 
00009 /* $Id: TOMS404.C,v 1.2.2.10 2008/07/17 17:07:08 garloff Exp $ */
00010 
00011 #if defined(__GNUC__) && __GNUC__ == 2 && defined(__cplusplus)
00012 # undef __cplusplus     /* workaround for buggy gcc-2.96 */
00013 #endif
00014 
00015 #include "prototypes.h"
00016 #include "prototypes2.h"
00017 
00018 #ifdef __cplusplus
00019 extern "C" {
00020 #endif  
00021 #include "f2c.h"
00022 /* Table of constant values */
00023 
00024 static integer c__2 = 2;
00025 
00026 /*     ALGORITHM 404 COLLECTED ALGORITHMS FROM ACM. */
00027 /*     ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 01, */
00028 /*     P. 048. */
00029 /* Complex */ /*VOID*/
00030 void cgamma_(complex *ret_val, complex *z__)
00031 {
00032     /* Format strings */
00033     static char fmt_900[] = "(1x,2e14.7,10x,\002ARGUMENT OF GAMMA FUNCTION I\
00034 S TOO CLOSE TO A POLE\002)";
00035     static char fmt_910[] = "(\002 ERROR - STIRLING*S SERIES HAS NOT CONVERG\
00036 ED\002/14x,\002Z = \002,2e14.7)";
00037 
00038     /* System generated locals */
00039     integer i__1;
00040     real r__1;
00041     complex q__1, q__2, q__3, q__4, q__5, q__6;
00042 
00043     /* Builtin functions: declared in prototypes.h */
00044 
00045     /* Local variables */
00046     static complex term;
00047     static integer iout;
00048     static complex a;
00049     static real c__[12];
00050     static integer i__, j, m;
00051     static complex t;
00052     static real x, y, xdist;
00053     static complex pi, zm, tt;
00054     static logical reflek;
00055     static complex den;
00056     static real tol;
00057     static complex sum;
00058 
00059     /* Fortran I/O blocks */
00060     static cilist io___9 = { 0, 0, 0, fmt_900, 0 };
00061     static cilist io___18 = { 0, 0, 0, fmt_910, 0 };
00062 
00063 
00064 /* SET IOUT FOR PROPER OUTPUT CHANNEL OF COMPUTER SYSTEM FOR */
00065 /* ERROR MESSAGES */
00066     iout = 3;
00067     pi.r = (float)3.14159265359, pi.i = (float)0.;
00068     x = z__->r;
00069     y = r_imag(z__);
00070 /* TOL = LIMIT OF PRECISION OF COMPUTER SYSTEM IN SINGLE PRECISI */
00071     tol = (float)1e-7;
00072     reflek = TRUE_;
00073 /* DETERMINE WHETHER Z IS TOO CLOSE TO A POLE */
00074 /* CHECK WHETHER TOO CLOSE TO ORIGIN */
00075     if (x >= tol) {
00076         goto L20;
00077     }
00078 /* FIND THE NEAREST POLE AND COMPUTE DISTANCE TO IT */
00079     xdist = x - (integer) (x - (float).5);
00080     q__1.r = xdist, q__1.i = y;
00081     zm.r = q__1.r, zm.i = q__1.i;
00082     if (c_abs(&zm) >= tol) {
00083         goto L10;
00084     }
00085 /* IF Z IS TOO CLOSE TO A POLE, PRINT ERROR MESSAGE AND RETURN */
00086 /* WITH CGAMMA = (1.E7,0.0E0) */
00087     io___9.ciunit = iout;
00088     s_wsfe(&io___9);
00089     do_fio(&c__2, (char *)&(*z__), (ftnlen)sizeof(real));
00090     e_wsfe();
00091      ret_val->r = (float)1e7,  ret_val->i = (float)0.;
00092     return ;
00093 /* FOR REAL(Z) NEGATIVE EMPLOY THE REFLECTION FORMULA */
00094 /* GAMMA(Z) = PI/(SIN(PI*Z)*GAMMA(1-Z)) */
00095 /* AND COMPUTE GAMMA(1-Z).  NOTE REFLEK IS A TAG TO INDICATE THA */
00096 /* THIS RELATION MUST BE USED LATER. */
00097 L10:
00098     if (x >= (float)0.) {
00099         goto L20;
00100     }
00101     reflek = FALSE_;
00102     q__1.r = (float)1. - z__->r, q__1.i = (float)0. - z__->i;
00103     z__->r = q__1.r, z__->i = q__1.i;
00104     x = (float)1. - x;
00105     y = -y;
00106 /* IF Z IS NOT TOO CLOSE TO A POLE, MAKE REAL(Z)>10 AND ARG(Z)<P */
00107 L20:
00108     m = 0;
00109 L40:
00110     if (x >= (float)10.) {
00111         goto L50;
00112     }
00113     x += (float)1.;
00114     ++m;
00115     goto L40;
00116 L50:
00117     if (dabs(y) < x) {
00118         goto L60;
00119     }
00120     x += (float)1.;
00121     ++m;
00122     goto L50;
00123 L60:
00124     q__1.r = x, q__1.i = y;
00125     t.r = q__1.r, t.i = q__1.i;
00126     q__1.r = t.r * t.r - t.i * t.i, q__1.i = t.r * t.i + t.i * t.r;
00127     tt.r = q__1.r, tt.i = q__1.i;
00128     den.r = t.r, den.i = t.i;
00129 /* COEFFICIENTS IN STIRLING*S APPROXIMATION FOR LN(GAMMA(T)) */
00130     c__[0] = (float).08333333333333333;
00131     c__[1] = (float)-.002777777777777778;
00132     c__[2] = (float)7.936507936507937e-4;
00133     c__[3] = (float)-5.952380952380953e-4;
00134     c__[4] = (float)8.417508417508417e-4;
00135     c__[5] = (float)-.001917526917526918;
00136     c__[6] = (float).00641025641025641;
00137     c__[7] = (float)-.02955065359477124;
00138     c__[8] = (float).1796443723688306;
00139     c__[9] = (float)-1.392432216905901;
00140     c__[10] = (float)13.40286404416839;
00141     q__4.r = t.r - (float).5, q__4.i = t.i + (float)0.;
00142     c_log(&q__5, &t);
00143     q__3.r = q__4.r * q__5.r - q__4.i * q__5.i, q__3.i = q__4.r * q__5.i + 
00144             q__4.i * q__5.r;
00145     q__2.r = q__3.r - t.r, q__2.i = q__3.i - t.i;
00146     r__1 = log((float)6.28318) * (float).5;
00147     q__6.r = r__1, q__6.i = (float)0.;
00148     q__1.r = q__2.r + q__6.r, q__1.i = q__2.i + q__6.i;
00149     sum.r = q__1.r, sum.i = q__1.i;
00150     j = 1;
00151 L70:
00152     i__1 = j - 1;
00153     q__2.r = c__[i__1], q__2.i = (float)0.;
00154     c_div(&q__1, &q__2, &den);
00155     term.r = q__1.r, term.i = q__1.i;
00156 /* TEST REAL AND IMAGINARY PARTS OF LN(GAMMA(Z)) SEPARATELY FOR */
00157 /* CONVERGENCE.  IF Z IS REAL SKIP IMAGINARY PART OF CHECK. */
00158     if ((r__1 = term.r / sum.r, dabs(r__1)) >= tol) {
00159         goto L80;
00160     }
00161     if (y == (float)0.) {
00162         goto L100;
00163     }
00164     if ((r__1 = r_imag(&term) / r_imag(&sum), dabs(r__1)) < tol) {
00165         goto L100;
00166     }
00167 L80:
00168     q__1.r = sum.r + term.r, q__1.i = sum.i + term.i;
00169     sum.r = q__1.r, sum.i = q__1.i;
00170     ++j;
00171     q__1.r = den.r * tt.r - den.i * tt.i, q__1.i = den.r * tt.i + den.i * 
00172             tt.r;
00173     den.r = q__1.r, den.i = q__1.i;
00174 /* TEST FOR NONCONVERGENCE */
00175     if (j == 12) {
00176         goto L90;
00177     }
00178     goto L70;
00179 /* STIRLING*S SERIES DID NOT CONVERGE.  PRINT ERROR MESSAGE AND */
00180 /* PROCEDE. */
00181 L90:
00182     io___18.ciunit = iout;
00183     s_wsfe(&io___18);
00184     do_fio(&c__2, (char *)&(*z__), (ftnlen)sizeof(real));
00185     e_wsfe();
00186 /* RECURSION RELATION USED TO OBTAIN LN(GAMMA(Z)) */
00187 /* LN(GAMMA(Z)) = LN(GAMMA(Z+M)/(Z*(Z+1)*...*(Z+M-1))) */
00188 /* = LN(GAMMA(Z+M)-LN(Z)-LN(Z+1)-...-LN(Z+M */
00189 L100:
00190     if (m == 0) {
00191         goto L120;
00192     }
00193     i__1 = m;
00194     for (i__ = 1; i__ <= i__1; ++i__) {
00195         r__1 = i__ * (float)1. - (float)1.;
00196         q__1.r = r__1, q__1.i = (float)0.;
00197         a.r = q__1.r, a.i = q__1.i;
00198 /* L110: */
00199         q__3.r = z__->r + a.r, q__3.i = z__->i + a.i;
00200         c_log(&q__2, &q__3);
00201         q__1.r = sum.r - q__2.r, q__1.i = sum.i - q__2.i;
00202         sum.r = q__1.r, sum.i = q__1.i;
00203     }
00204 /* CHECK TO SEE IF REFLECTION FORMULA SHOULD BE USED */
00205 L120:
00206     if (reflek) {
00207         goto L130;
00208     }
00209     q__5.r = pi.r * z__->r - pi.i * z__->i, q__5.i = pi.r * z__->i + pi.i * 
00210             z__->r;
00211     c_sin(&q__4, &q__5);
00212     c_div(&q__3, &pi, &q__4);
00213     c_log(&q__2, &q__3);
00214     q__1.r = q__2.r - sum.r, q__1.i = q__2.i - sum.i;
00215     sum.r = q__1.r, sum.i = q__1.i;
00216     q__1.r = (float)1. - z__->r, q__1.i = (float)0. - z__->i;
00217     z__->r = q__1.r, z__->i = q__1.i;
00218 L130:
00219     c_exp(&q__1, &sum);
00220      ret_val->r = q__1.r,  ret_val->i = q__1.i;
00221     return ;
00222 } /* cgamma_ */
00223 
00224 #ifdef __cplusplus
00225 }
00226 #endif
00227