1      SUBROUTINE  ZGAM(CARG,CANS,ERREST,MODE)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5c>> 1996-03-30 ZGAM Krogh Added external statement.
6C>> 1995-11-20 ZGAM Krogh Set up so M77CON converts between "Z" and "C".
7C>> 1994-08-17 CLL Add tests on BIGINT to allow easier conversion to C.
8C>> 1994-05-25 ZGAM WVS generate COEF using PARAMETER
9C>> 1994-04-20 ZGAM CLL Make DP and SP versions similar.
10C>> 1993-04-13 ZGAM CLL Edit for conversion to C.
11C>> 1992-04-20 ZGAM CLL Edited comments.
12C>> 1991-11-11 ZGAM CLL Made [Z/C]GAM from CDLGAM
13C>> 1991-01-16 CDLGAM Lawson  Removing use of subr D2MACH.
14C>> 1985-08-02 CDLGAM Lawson  Initial code.
15C
16C *** COMPLEX GAMMA AND LOGGAMMA FUNCTIONS WITH ERROR ESTIMATE
17C
18C     -----------------------------------------------------------------
19C     SUBROUTINE ARGUMENTS
20C     --------------------
21C     CARG()   A complex argument, given as an array of 2 floating-point
22c              elements consisting of the real component
23c              followed by the imaginary component.
24c
25c     CANS()   The complex answer, stored as an array of 2
26c              floating-point numbers, representing the real and
27c              imaginary parts.
28c
29c     ERREST   On output ERREST gives an estimate of the absolute
30c              (for LOGGAMMA) or the relative (for GAMMA) error
31c              of the answer.
32c
33c     MODE     Selects function to be computed.  set it to 0 for
34c              LOGGAMMA, and 1 for GAMMA.
35C     -----------------------------------------------------------------
36C     MACHINE DEPENDANT PARAMETERS
37C     If the fraction part of a floating point number
38C     contains T digits using base B then
39C         EPS3  = B ** (-T)
40C         EPS4  = B ** (-T+1)
41C         OMEGA = overflow limit
42C         DESET = 5.0 on a binary machine
43C               = 2.0 on a base 16 machine
44C     -----------------------------------------------------------------
45C     REFERENCE: H.KUKI, Comm.ACM, Vol.15, (1972),
46C     pp.262-267, 271-272.  Subroutine name was CDLGAM.
47C     Code developed for UNIVAC 1108 by E.W.NG, JPL, 1969.
48C     Modified for FORTRAN 77 portability by C.L.LAWSON &
49C     S.CHAN, JPL, 1983.
50c     -----------------------------------------------------------------
51c--Z replaces "?": ?GAM
52c--D (type)replaces "?": ?ERM1, ?ERV1
53c     Also uses I1MACH, and D1MACH
54c     -----------------------------------------------------------------
55      external D1MACH, I1MACH
56      double precision D1MACH
57      double precision A,AL1,AL2,AL2P,B,BIGINT
58      double precision CARG(2),CANS(2),COEF(7), CUT1
59      double precision DE0,DE1,DELTA,DESET,DN
60      double precision ELIMIT, EPS3, EPS4,ERREST,F0,F1,G0,G1
61      double precision H,H1,H2,HL2P,OMEGA, ONE
62      double precision PI, REPS3, T1,T2,TWOPI
63      double precision U,U1,U2,UU1,UU2,UUU1,UUU2,V1,V2,VV1,VV2
64      double precision W1,W2,W3,Y1,Z1,Z2, ZERO, ZZ1
65      integer I1MACH, ITEMP, J, K, LF1, LF2, LF3, MODE, N
66      logical FIRST
67C
68      SAVE FIRST,BIGINT,COEF,OMEGA,EPS4,EPS3,REPS3,CUT1,DESET,ELIMIT
69C
70      parameter(ONE = 1.0D0, ZERO = 0.0D0)
71      parameter(F0 = 840.07385296052619D0, F1 = 20.001230821894200D0)
72      parameter(G0 = 1680.1477059210524D0, G1 = 180.01477047052042D0)
73      parameter(PI = 3.14159265358979323846D0)
74      parameter(TWOPI = 6.283185307179586476925D0)
75      parameter(HL2P = 0.918938533204672742D0)
76      parameter(AL2P = 1.83787706640934548D0)
77
78c              COEF(8-i) = bernoulli(2i)/(2i*(2i-1)).
79      double precision C1, C2, C3, C4, C5, C6, C7
80      parameter (C1 = +1.0e0/156.0e0)
81      parameter (C2 = -691.0e0/360360.0e0)
82      parameter (C3 = +1.0e0/1188.0e0)
83      parameter (C4 = -1.0e0/1680.0e0)
84      parameter (C5 = +1.0e0/1260.0e0)
85      parameter (C6 = -1.0e0/360.0e0)
86      parameter (C7 = +1.0e0/12.0e0)
87      DATA COEF / C1, C2, C3, C4, C5, C6, C7 /
88c     DATA COEF   /+0.641025641025641026E-2,
89c    *             -0.191752691752691753E-2,
90c    *             +0.841750841750841751E-3,
91c    *             -0.595238095238095238E-3,
92c    *             +0.793650793650793651E-3,
93c    *             -0.277777777777777778E-2,
94c    *             +0.833333333333333333E-1/
95      DATA FIRST/.TRUE./
96c     ------------------------------------------------------------------
97      IF (FIRST) THEN
98        FIRST = .FALSE.
99        OMEGA = D1MACH(2)
100        EPS3  = D1MACH(3)
101        EPS4  = D1MACH(4)
102        REPS3 = ONE / EPS3
103        ELIMIT = log(OMEGA)
104        CUT1  = log(EPS3)
105        BIGINT = I1MACH(9) - 2
106          IF (I1MACH(10) .EQ. 2) THEN
107            DESET = 5.D0
108          ELSE
109            DESET = 2.D0
110          ENDIF
111      ENDIF
112      DE0 = DESET
113      DE1 = ZERO
114      Z1 = CARG(1)
115      Z2 = CARG(2)
116C
117C *** Setting DELTA = estimate of uncertainty level of
118C                     argument data.
119C
120      DELTA = EPS4 * (abs(Z1) + abs(Z2))
121      IF (DELTA .EQ. ZERO) DELTA = EPS4
122C
123C *** FORCE SIGN OF IMAGINARY PART OF ARG TO NON-NEGATIVE
124C
125      LF1 = 0
126      IF (Z2 .lt. ZERO) then
127         LF1 = 1
128         Z2 = -Z2
129      endif
130      LF2 = 0
131      IF (Z1 .GE. ZERO) GO TO 100
132C
133C *** CASE WHEN REAL PART OF ARG IS NEGATIVE
134C
135      LF2 = 1
136      LF1 = LF1-1
137      T1 = AL2P - PI*Z2
138      T2 = PI*(0.5D0 - Z1)
139      U  = -TWOPI*Z2
140      IF (U .lt. -0.1054D0) then
141         A  = ZERO
142C
143C ***    IF EXP(U) .LE. EPS3, IGNORE IT TO SAVE TIME AND TO AVOID
144C        IRRELEVANT UNDERFLOW
145C
146         IF (U .gt. CUT1) then
147            A  = exp(U)
148         endif
149         H1 = ONE - A
150      else
151         U2 = U * U
152         A  = -U*(F1*U2 + F0)
153         H1 = (A + A)/((U2 + G1)*U2 + G0 + A)
154         A  = ONE - H1
155      endif
156c
157c              Here Z1 is negative.
158c
159      if(Z1 .lt. -BIGINT) then
160         CALL DERM1('ZGAM',3,0,'Require CARG(1) .ge. -BIGINT',
161     *              'CARG(1)', Z1, ',')
162         CALL DERV1('-BIGINT',-BIGINT,'.')
163         go to 700
164      endif
165c
166c           Truncate to integer: ITEMP
167c
168      ITEMP = Z1 - 0.5d0
169      B  = Z1 - ITEMP
170      H2 = A*sin(TWOPI*B)
171      B  = sin(PI*B)
172      H1 = H1 + (B+B)*B*A
173      H  = abs(H2) + H1 - TWOPI*A*DELTA
174      IF (H .LE. ZERO) GO TO 500
175      DE0 = DE0 + abs(T1) + T2
176      DE1 = PI + TWOPI*A/H
177      Z1 = ONE - Z1
178C
179C *** CASE WHEN NEITHER REAL PART NOR IMAGINARY PART OF ARG IS
180C     NEGATIVE.  DEFINE THRESHOLD CURVE TO BE THE BROKEN LINES
181C     CONNECTING POINTS 10F0*I, 10F4.142*I, 0.1F14.042*I,AND
182C     0.1FOMEGA*I
183C
184  100 LF3 = 0
185      Y1 = Z1 - 0.5D0
186      W1 = ZERO
187      W2 = ZERO
188      K  = 0
189      B  = max(0.1D0, min(10.0D0, 14.142D0-Z2)) - Z1
190      IF (B .LE. ZERO) GO TO 200
191C
192C *** CASE WHEN REAL PART OF ARG IS BETWEEN 0 AND THRESHOLD
193C
194      LF3 = 1
195      ZZ1 = Z1
196      N  = B + ONE
197      DN = N
198      Z1 = Z1 + DN
199      A  = Z1*Z1 + Z2*Z2
200      V1 = Z1/A
201      V2 = -Z2/A
202C
203C *** INITIALIZE U1+U2*I AS THE RIGHTMOST FACTOR 1-1/(Z+N)
204C
205      U1 = ONE - V1
206      U2 = -V2
207      K  = 6.0D0 - Z2*0.6D0 - ZZ1
208      IF (K .gt. 0) then
209C
210C *** FORWARD ASSEMBLY OF FACTORS (Z+J-1)/(Z+N)
211C
212         N  = N - K
213         UU1 = (ZZ1*Z1 + Z2*Z2) / A
214         UU2 = DN*Z2/A
215         VV1 = ZERO
216         VV2 = ZERO
217         DO 110  J = 1,K
218            B  = U1*(UU1+VV1) - U2*(UU2+VV2)
219            U2 = U1*(UU2+VV2) + U2*(UU1+VV1)
220            U1 = B
221            VV1 = VV1 + V1
222            VV2 = VV2 + V2
223  110    continue
224      endif
225      IF (N .GE. 2) then
226C
227C *** BACKWARD ASSEMBLY OF FACTORS 1-J/(Z+N)
228C
229         VV1 = V1
230         VV2 = V2
231         DO 130  J = 2,N
232            VV1 = VV1 + V1
233            VV2 = VV2 + V2
234            B  = U1*(ONE - VV1) + U2*VV2
235            U2 = -U1*VV2 + U2*(ONE - VV1)
236            U1 = B
237  130   continue
238      endif
239      U  = U1*U1 + U2*U2
240      IF (U .EQ. ZERO) GO TO 500
241      IF (MODE .NE. 0) then
242         IF (K .LE. 0) GO TO 200
243      endif
244      AL1 = log(U)*0.5D0
245      IF (MODE .EQ. 0) then
246         W1 = AL1
247         W2 = atan2(U2,U1)
248         IF (W2 .LT. ZERO)  W2 = W2 + TWOPI
249         IF (K .LE. 0) GO TO 200
250      endif
251      A  = ZZ1 + Z2 - DELTA
252      IF (A .LE. ZERO) GO TO 500
253      DE0 = DE0 - AL1
254      DE1 = DE1 + 2.0D0 + ONE/A
255C
256C *** CASE WHEN REAL PART OF ARG IS GREATER THAN THRESHOLD
257C
258  200 A  = Z1*Z1 + Z2*Z2
259      AL1 = log(A)*0.5D0
260      AL2 = atan2(Z2,Z1)
261      V1 = Y1*AL1 - Z2*AL2
262      V2 = Y1*AL2 + Z2*AL1
263C
264C *** Evaluate asymptotic terms. Ignore this term,
265C     if ABS(ARG)**2 .GT. REPS3, to save time and
266C     to avoid irrelevant underflow.
267C
268      VV1 = ZERO
269      VV2 = ZERO
270      IF (A .GT. REPS3) GO TO 220
271      UU1 = Z1/A
272      UU2 = -Z2/A
273      UUU1 = UU1*UU1 - UU2*UU2
274      UUU2 = UU1*UU2*2.0D0
275      VV1 = COEF(1)
276      DO 210  J = 2,7
277        B  = VV1*UUU1 - VV2*UUU2
278        VV2 = VV1*UUU2 + VV2*UUU1
279  210   VV1 = B + COEF(J)
280      B  = VV1*UU1 - VV2*UU2
281      VV2 = VV1*UU2 + VV2*UU1
282      VV1 = B
283  220 W1 = (((VV1 + HL2P) - W1) - Z1) + V1
284      W2 = ((VV2 - W2) - Z2) + V2
285      DE0 = DE0 + abs(V1) + abs(V2)
286      IF (K .LE. 0)  DE1 = DE1 + AL1
287C FINAL ASSEMBLY
288      IF (LF2 .EQ. 0) then
289         IF (MODE .NE. 0) then
290            IF (W1 .GT. ELIMIT) GO TO 550
291            A  = exp(W1)
292            W1 = A*cos(W2)
293            W2 = A*sin(W2)
294            IF (LF3 .NE. 0) then
295               B  = (W1*U1 + W2*U2) / U
296               W2 = (W2*U1 - W1*U2) / U
297               W1 = B
298            endif
299         endif
300         GO TO 400
301      endif
302      H  = H1*H1 + H2*H2
303      IF (H .EQ. ZERO) GO TO 500
304      IF (MODE .EQ. 0 .or. H .le. 1.0D-2) then
305         A  = log(H)*0.5D0
306         IF (H .LE. 1.0D-2)  DE0 = DE0 - A
307         IF (MODE .EQ. 0) then
308            W1 = (T1 - A) - W1
309            W2 = (T2 - atan2(H2,H1)) - W2
310            GO TO 400
311         endif
312      endif
313c                       Here we have MODE .ne. 0 and LF2 .ne. 0.
314      T1 = T1 - W1
315      T2 = T2 - W2
316      IF (T1 .GT. ELIMIT) GO TO 550
317      A  = exp(T1)
318      T1 = A*cos(T2)
319      T2 = A*sin(T2)
320      W1 = (T1*H1 + T2*H2)/H
321      W2 = (T2*H1 - T1*H2)/H
322      IF (LF3 .NE. 0) then
323         B  = W1*U1 - W2*U2
324         W2 = W1*U2 + W2*U1
325         W1 = B
326      endif
327  400 continue
328      IF (LF1 .NE. 0)  W2 = -W2
329C
330C *** TRUNCATION ERREST OF STIRLINGS FORMULA IS UP TO EPS3.
331C
332      DE1 = DE0*EPS4 + EPS3 + DE1*DELTA
333
334c                                      Normal termination.
335c
336c     The imaginary part of the log of a complex number is nonunique
337c     to within multiples of 2*Pi.  We prefer a result for loggamma
338c     having its imaginary part .gt. -Pi and .le. +Pi.  The result at
339c     this point is usually in this range.  If not we will move it
340c     into this range.  -- CLL 11/11/91
341c
342      if(MODE .eq. 0) then
343         if(W2 .le. -PI .or. W2 .gt. PI) then
344            W3 = abs(W2)
345            T1 = W3/PI + ONE
346            if(abs(T1) .gt. BIGINT) then
347               CALL DERM1('ZGAM',4,0,'Argument out of range.',
348     *                    'CARG(1)',CARG(1),',')
349               CALL DERV1('CARG(2)',CARG(2),'.')
350               go to 700
351            endif
352            T2 = int(T1) / 2
353            W3 = W3 - T2 * TWOPI
354            if(W2 .ge. 0.0d0) then
355               W2 = W3
356            else
357               W2 = -W3
358            endif
359            if(W2 .le. -PI) then
360               W2 = W2 + TWOPI
361            elseif(W2 .gt. PI) then
362               W2 = W2 - TWOPI
363            endif
364         endif
365      endif
366      CANS(1) = W1
367      CANS(2) = W2
368      ERREST = DE1
369      return
370c                                           Error termination.
371C
372C *** CASE WHEN ARGUMENT IS TOO CLOSE TO A SINGULARITY
373C
374  500 CONTINUE
375      CALL DERM1('ZGAM',1,0,'Z TOO CLOSE TO A SINGULARITY',
376     *           'Z(1)',CARG(1),',')
377      CALL DERV1('Z(2)',CARG(2),'.')
378      GO TO 700
379C
380  550 CONTINUE
381      CALL DERM1('ZGAM',2,0,'ARG TOO LARGE. EXP FUNCTION OVERFLOW',
382     *           'Z(1)',CARG(1),',')
383      CALL DERV1('Z(2)',CARG(2),'.')
384  700 continue
385      CANS(1) = OMEGA
386      CANS(2) = OMEGA
387      ERREST = OMEGA
388      return
389      end
390