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