/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:49 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "zgam.h" #include #include #include /* PARAMETER translations */ #define AL2P 1.83787706640934548e0 #define C1 (1.0e0/156.0e0) #define C2 (-691.0e0/360360.0e0) #define C3 (1.0e0/1188.0e0) #define C4 (-1.0e0/1680.0e0) #define C5 (1.0e0/1260.0e0) #define C6 (-1.0e0/360.0e0) #define C7 (1.0e0/12.0e0) #define F0 840.07385296052619e0 #define F1 20.001230821894200e0 #define G0 1680.1477059210524e0 #define G1 180.01477047052042e0 #define HL2P 0.918938533204672742e0 #define ONE 1.0e0 #define PI 3.14159265358979323846e0 #define TWOPI 6.283185307179586476925e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ zgam( double carg[], double cans[], double *errest, long mode) { long int _l0, itemp, j, k, lf1, lf2, lf3, n; double a, al1, al2, b, de0, de1, delta, dn, h, h1, h2, t1, t2, u, u1, u2, uu1, uu2, uuu1, uuu2, v1, v2, vv1, vv2, w1, w2, w3, y1, z1, z2, zz1; static double bigint, cut1, deset, elimit, eps3, eps4, omega, reps3; static double coef[7]={C1,C2,C3,C4,C5,C6,C7}; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Cans = &cans[0] - 1; double *const Carg = &carg[0] - 1; double *const Coef = &coef[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-03-30 ZGAM Krogh Added external statement. *>> 1995-11-20 ZGAM Krogh Set up so M77CON converts between "Z" and "C". *>> 1994-08-17 CLL Add tests on BIGINT to allow easier conversion to C. *>> 1994-05-25 ZGAM WVS generate COEF using PARAMETER *>> 1994-04-20 ZGAM CLL Make DP and SP versions similar. *>> 1993-04-13 ZGAM CLL Edit for conversion to C. *>> 1992-04-20 ZGAM CLL Edited comments. *>> 1991-11-11 ZGAM CLL Made [Z/C]GAM from CDLGAM *>> 1991-01-16 CDLGAM Lawson Removing use of subr D2MACH. *>> 1985-08-02 CDLGAM Lawson Initial code. * * *** COMPLEX GAMMA AND LOGGAMMA FUNCTIONS WITH ERROR ESTIMATE * * ----------------------------------------------------------------- * SUBROUTINE ARGUMENTS * -------------------- * CARG() A complex argument, given as an array of 2 floating-point * elements consisting of the real component * followed by the imaginary component. * * CANS() The complex answer, stored as an array of 2 * floating-point numbers, representing the real and * imaginary parts. * * ERREST On output ERREST gives an estimate of the absolute * (for LOGGAMMA) or the relative (for GAMMA) error * of the answer. * * MODE Selects function to be computed. set it to 0 for * LOGGAMMA, and 1 for GAMMA. * ----------------------------------------------------------------- * MACHINE DEPENDANT PARAMETERS * If the fraction part of a floating point number * contains T digits using base B then * EPS3 = B ** (-T) * EPS4 = B ** (-T+1) * OMEGA = overflow limit * DESET = 5.0 on a binary machine * = 2.0 on a base 16 machine * ----------------------------------------------------------------- * REFERENCE: H.KUKI, Comm.ACM, Vol.15, (1972), * pp.262-267, 271-272. Subroutine name was CDLGAM. * Code developed for UNIVAC 1108 by E.W.NG, JPL, 1969. * Modified for FORTRAN 77 portability by C.L.LAWSON & * S.CHAN, JPL, 1983. * ----------------------------------------------------------------- *--Z replaces "?": ?GAM *--D (type)replaces "?": ?ERM1, ?ERV1 * Also uses I1MACH, and D1MACH * ----------------------------------------------------------------- */ /* COEF(8-i) = bernoulli(2i)/(2i*(2i-1)). */ /* DATA COEF /+0.641025641025641026E-2, * * -0.191752691752691753E-2, * * +0.841750841750841751E-3, * * -0.595238095238095238E-3, * * +0.793650793650793651E-3, * * -0.277777777777777778E-2, * * +0.833333333333333333E-1/ */ /* ------------------------------------------------------------------ */ if (first) { first = FALSE; omega = DBL_MAX; eps3 = DBL_EPSILON/FLT_RADIX; eps4 = DBL_EPSILON; reps3 = ONE/eps3; elimit = log( omega ); cut1 = log( eps3 ); bigint = LONG_MAX - 2; if (FLT_RADIX == 2) { deset = 5.e0; } else { deset = 2.e0; } } de0 = deset; de1 = ZERO; z1 = Carg[1]; z2 = Carg[2]; /* *** Setting DELTA = estimate of uncertainty level of * argument data. * */ delta = eps4*(fabs( z1 ) + fabs( z2 )); if (delta == ZERO) delta = eps4; /* *** FORCE SIGN OF IMAGINARY PART OF ARG TO NON-NEGATIVE * */ lf1 = 0; if (z2 < ZERO) { lf1 = 1; z2 = -z2; } lf2 = 0; if (z1 >= ZERO) goto L_100; /* *** CASE WHEN REAL PART OF ARG IS NEGATIVE * */ lf2 = 1; lf1 -= 1; t1 = AL2P - PI*z2; t2 = PI*(0.5e0 - z1); u = -TWOPI*z2; if (u < -0.1054e0) { a = ZERO; /* *** IF EXP(U) .LE. EPS3, IGNORE IT TO SAVE TIME AND TO AVOID * IRRELEVANT UNDERFLOW * */ if (u > cut1) { a = exp( u ); } h1 = ONE - a; } else { u2 = u*u; a = -u*(F1*u2 + F0); h1 = (a + a)/((u2 + G1)*u2 + G0 + a); a = ONE - h1; } /* Here Z1 is negative. * */ if (z1 < -bigint) { derm1( "ZGAM", 3, 0, "Require CARG(1) .ge. -BIGINT", "CARG(1)" , z1, ',' ); derv1( "-BIGINT", -bigint, '.' ); goto L_700; } /* Truncate to integer: ITEMP * */ itemp = z1 - 0.5e0; b = z1 - itemp; h2 = a*sin( TWOPI*b ); b = sin( PI*b ); h1 += (b + b)*b*a; h = fabs( h2 ) + h1 - TWOPI*a*delta; if (h <= ZERO) goto L_500; de0 += fabs( t1 ) + t2; de1 = PI + TWOPI*a/h; z1 = ONE - z1; /* *** CASE WHEN NEITHER REAL PART NOR IMAGINARY PART OF ARG IS * NEGATIVE. DEFINE THRESHOLD CURVE TO BE THE BROKEN LINES * CONNECTING POINTS 10F0*I, 10F4.142*I, 0.1F14.042*I,AND * 0.1FOMEGA*I * */ L_100: lf3 = 0; y1 = z1 - 0.5e0; w1 = ZERO; w2 = ZERO; k = 0; b = fmax( 0.1e0, fmin( 10.0e0, 14.142e0 - z2 ) ) - z1; if (b <= ZERO) goto L_200; /* *** CASE WHEN REAL PART OF ARG IS BETWEEN 0 AND THRESHOLD * */ lf3 = 1; zz1 = z1; n = b + ONE; dn = n; z1 += dn; a = z1*z1 + z2*z2; v1 = z1/a; v2 = -z2/a; /* *** INITIALIZE U1+U2*I AS THE RIGHTMOST FACTOR 1-1/(Z+N) * */ u1 = ONE - v1; u2 = -v2; k = 6.0e0 - z2*0.6e0 - zz1; if (k > 0) { /* *** FORWARD ASSEMBLY OF FACTORS (Z+J-1)/(Z+N) * */ n -= k; uu1 = (zz1*z1 + z2*z2)/a; uu2 = dn*z2/a; vv1 = ZERO; vv2 = ZERO; for (j = 1; j <= k; j++) { b = u1*(uu1 + vv1) - u2*(uu2 + vv2); u2 = u1*(uu2 + vv2) + u2*(uu1 + vv1); u1 = b; vv1 += v1; vv2 += v2; } } if (n >= 2) { /* *** BACKWARD ASSEMBLY OF FACTORS 1-J/(Z+N) * */ vv1 = v1; vv2 = v2; for (j = 2; j <= n; j++) { vv1 += v1; vv2 += v2; b = u1*(ONE - vv1) + u2*vv2; u2 = -u1*vv2 + u2*(ONE - vv1); u1 = b; } } u = u1*u1 + u2*u2; if (u == ZERO) goto L_500; if (mode != 0) { if (k <= 0) goto L_200; } al1 = log( u )*0.5e0; if (mode == 0) { w1 = al1; w2 = atan2( u2, u1 ); if (w2 < ZERO) w2 += TWOPI; if (k <= 0) goto L_200; } a = zz1 + z2 - delta; if (a <= ZERO) goto L_500; de0 -= al1; de1 += 2.0e0 + ONE/a; /* *** CASE WHEN REAL PART OF ARG IS GREATER THAN THRESHOLD * */ L_200: a = z1*z1 + z2*z2; al1 = log( a )*0.5e0; al2 = atan2( z2, z1 ); v1 = y1*al1 - z2*al2; v2 = y1*al2 + z2*al1; /* *** Evaluate asymptotic terms. Ignore this term, * if ABS(ARG)**2 .GT. REPS3, to save time and * to avoid irrelevant underflow. * */ vv1 = ZERO; vv2 = ZERO; if (a > reps3) goto L_220; uu1 = z1/a; uu2 = -z2/a; uuu1 = uu1*uu1 - uu2*uu2; uuu2 = uu1*uu2*2.0e0; vv1 = Coef[1]; for (j = 2; j <= 7; j++) { b = vv1*uuu1 - vv2*uuu2; vv2 = vv1*uuu2 + vv2*uuu1; vv1 = b + Coef[j]; } b = vv1*uu1 - vv2*uu2; vv2 = vv1*uu2 + vv2*uu1; vv1 = b; L_220: w1 = (((vv1 + HL2P) - w1) - z1) + v1; w2 = ((vv2 - w2) - z2) + v2; de0 += fabs( v1 ) + fabs( v2 ); if (k <= 0) de1 += al1; /* FINAL ASSEMBLY */ if (lf2 == 0) { if (mode != 0) { if (w1 > elimit) goto L_550; a = exp( w1 ); w1 = a*cos( w2 ); w2 = a*sin( w2 ); if (lf3 != 0) { b = (w1*u1 + w2*u2)/u; w2 = (w2*u1 - w1*u2)/u; w1 = b; } } goto L_400; } h = h1*h1 + h2*h2; if (h == ZERO) goto L_500; if (mode == 0 || h <= 1.0e-2) { a = log( h )*0.5e0; if (h <= 1.0e-2) de0 -= a; if (mode == 0) { w1 = (t1 - a) - w1; w2 = (t2 - atan2( h2, h1 )) - w2; goto L_400; } } /* Here we have MODE .ne. 0 and LF2 .ne. 0. */ t1 -= w1; t2 -= w2; if (t1 > elimit) goto L_550; a = exp( t1 ); t1 = a*cos( t2 ); t2 = a*sin( t2 ); w1 = (t1*h1 + t2*h2)/h; w2 = (t2*h1 - t1*h2)/h; if (lf3 != 0) { b = w1*u1 - w2*u2; w2 = w1*u2 + w2*u1; w1 = b; } L_400: ; if (lf1 != 0) w2 = -w2; /* *** TRUNCATION ERREST OF STIRLINGS FORMULA IS UP TO EPS3. * */ de1 = de0*eps4 + eps3 + de1*delta; /* Normal termination. * * The imaginary part of the log of a complex number is nonunique * to within multiples of 2*Pi. We prefer a result for loggamma * having its imaginary part .gt. -Pi and .le. +Pi. The result at * this point is usually in this range. If not we will move it * into this range. -- CLL 11/11/91 * */ if (mode == 0) { if (w2 <= -PI || w2 > PI) { w3 = fabs( w2 ); t1 = w3/PI + ONE; if (fabs( t1 ) > bigint) { derm1( "ZGAM", 4, 0, "Argument out of range.", "CARG(1)" , Carg[1], ',' ); derv1( "CARG(2)", Carg[2], '.' ); goto L_700; } t2 = (long)( t1 )/2; w3 += -t2*TWOPI; if (w2 >= 0.0e0) { w2 = w3; } else { w2 = -w3; } if (w2 <= -PI) { w2 += TWOPI; } else if (w2 > PI) { w2 -= TWOPI; } } } Cans[1] = w1; Cans[2] = w2; *errest = de1; return; /* Error termination. * * *** CASE WHEN ARGUMENT IS TOO CLOSE TO A SINGULARITY * */ L_500: ; derm1( "ZGAM", 1, 0, "Z TOO CLOSE TO A SINGULARITY", "Z(1)", Carg[1], ',' ); derv1( "Z(2)", Carg[2], '.' ); goto L_700; L_550: ; derm1( "ZGAM", 2, 0, "ARG TOO LARGE. EXP FUNCTION OVERFLOW", "Z(1)" , Carg[1], ',' ); derv1( "Z(2)", Carg[2], '.' ); L_700: ; Cans[1] = omega; Cans[2] = omega; *errest = omega; return; } /* end of function */