github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zunk1.f (about) 1 SUBROUTINE ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, 2 * ALIM) 3 C***BEGIN PROLOGUE ZUNK1 4 C***REFER TO ZBESK 5 C 6 C ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE 7 C RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE 8 C UNIFORM ASYMPTOTIC EXPANSION. 9 C MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. 10 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR 11 C 12 C***ROUTINES CALLED ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,ZABS 13 C***END PROLOGUE ZUNK1 14 C COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO, 15 C *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR 16 DOUBLE PRECISION ALIM, ANG, APHI, ASC, ASCLE, BRY, CKI, CKR, 17 * CONER, CRSC, CSCL, CSGNI, CSPNI, CSPNR, CSR, CSRR, CSSR, 18 * CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2M, C2R, ELIM, FMR, FN, 19 * FNF, FNU, PHIDI, PHIDR, PHII, PHIR, PI, RAST, RAZR, RS1, RZI, 20 * RZR, SGN, STI, STR, SUMDI, SUMDR, SUMI, SUMR, S1I, S1R, S2I, 21 * S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, 22 * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, ZABS 23 INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG, 24 * KK, KODE, MR, N, NW, NZ, INITD, IC, IPARD, J 25 DIMENSION BRY(3), INIT(2), YR(N), YI(N), SUMR(2), SUMI(2), 26 * ZETA1R(2), ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2), 27 * CWRKR(16,3), CWRKI(16,3), CSSR(3), CSRR(3), PHIR(2), PHII(2) 28 DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / 29 DATA PI / 3.14159265358979324D0 / 30 C 31 KDFLG = 1 32 NZ = 0 33 C----------------------------------------------------------------------- 34 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN 35 C THE UNDERFLOW LIMIT 36 C----------------------------------------------------------------------- 37 CSCL = 1.0D0/TOL 38 CRSC = TOL 39 CSSR(1) = CSCL 40 CSSR(2) = CONER 41 CSSR(3) = CRSC 42 CSRR(1) = CRSC 43 CSRR(2) = CONER 44 CSRR(3) = CSCL 45 BRY(1) = 1.0D+3*D1MACH(1)/TOL 46 BRY(2) = 1.0D0/BRY(1) 47 BRY(3) = D1MACH(2) 48 ZRR = ZR 49 ZRI = ZI 50 IF (ZR.GE.0.0D0) GO TO 10 51 ZRR = -ZR 52 ZRI = -ZI 53 10 CONTINUE 54 J = 2 55 DO 70 I=1,N 56 C----------------------------------------------------------------------- 57 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J 58 C----------------------------------------------------------------------- 59 J = 3 - J 60 FN = FNU + DBLE(FLOAT(I-1)) 61 INIT(J) = 0 62 CALL ZUNIK(ZRR, ZRI, FN, 2, 0, TOL, INIT(J), PHIR(J), PHII(J), 63 * ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), SUMR(J), SUMI(J), 64 * CWRKR(1,J), CWRKI(1,J)) 65 IF (KODE.EQ.1) GO TO 20 66 STR = ZRR + ZETA2R(J) 67 STI = ZRI + ZETA2I(J) 68 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 69 STR = STR*RAST*RAST 70 STI = -STI*RAST*RAST 71 S1R = ZETA1R(J) - STR 72 S1I = ZETA1I(J) - STI 73 GO TO 30 74 20 CONTINUE 75 S1R = ZETA1R(J) - ZETA2R(J) 76 S1I = ZETA1I(J) - ZETA2I(J) 77 30 CONTINUE 78 RS1 = S1R 79 C----------------------------------------------------------------------- 80 C TEST FOR UNDERFLOW AND OVERFLOW 81 C----------------------------------------------------------------------- 82 IF (DABS(RS1).GT.ELIM) GO TO 60 83 IF (KDFLG.EQ.1) KFLAG = 2 84 IF (DABS(RS1).LT.ALIM) GO TO 40 85 C----------------------------------------------------------------------- 86 C REFINE TEST AND SCALE 87 C----------------------------------------------------------------------- 88 APHI = ZABS(CMPLX(PHIR(J),PHII(J),kind=KIND(1.0D0))) 89 RS1 = RS1 + DLOG(APHI) 90 IF (DABS(RS1).GT.ELIM) GO TO 60 91 IF (KDFLG.EQ.1) KFLAG = 1 92 IF (RS1.LT.0.0D0) GO TO 40 93 IF (KDFLG.EQ.1) KFLAG = 3 94 40 CONTINUE 95 C----------------------------------------------------------------------- 96 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR 97 C EXPONENT EXTREMES 98 C----------------------------------------------------------------------- 99 S2R = PHIR(J)*SUMR(J) - PHII(J)*SUMI(J) 100 S2I = PHIR(J)*SUMI(J) + PHII(J)*SUMR(J) 101 STR = DEXP(S1R)*CSSR(KFLAG) 102 S1R = STR*DCOS(S1I) 103 S1I = STR*DSIN(S1I) 104 STR = S2R*S1R - S2I*S1I 105 S2I = S1R*S2I + S2R*S1I 106 S2R = STR 107 IF (KFLAG.NE.1) GO TO 50 108 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 109 IF (NW.NE.0) GO TO 60 110 50 CONTINUE 111 CYR(KDFLG) = S2R 112 CYI(KDFLG) = S2I 113 YR(I) = S2R*CSRR(KFLAG) 114 YI(I) = S2I*CSRR(KFLAG) 115 IF (KDFLG.EQ.2) GO TO 75 116 KDFLG = 2 117 GO TO 70 118 60 CONTINUE 119 IF (RS1.GT.0.0D0) GO TO 300 120 C----------------------------------------------------------------------- 121 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 122 C----------------------------------------------------------------------- 123 IF (ZR.LT.0.0D0) GO TO 300 124 KDFLG = 1 125 YR(I)=ZEROR 126 YI(I)=ZEROI 127 NZ=NZ+1 128 IF (I.EQ.1) GO TO 70 129 IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 70 130 YR(I-1)=ZEROR 131 YI(I-1)=ZEROI 132 NZ=NZ+1 133 70 CONTINUE 134 I = N 135 75 CONTINUE 136 RAZR = 1.0D0/ZABS(CMPLX(ZRR,ZRI,kind=KIND(1.0D0))) 137 STR = ZRR*RAZR 138 STI = -ZRI*RAZR 139 RZR = (STR+STR)*RAZR 140 RZI = (STI+STI)*RAZR 141 CKR = FN*RZR 142 CKI = FN*RZI 143 IB = I + 1 144 IF (N.LT.IB) GO TO 160 145 C----------------------------------------------------------------------- 146 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO 147 C ON UNDERFLOW. 148 C----------------------------------------------------------------------- 149 FN = FNU + DBLE(FLOAT(N-1)) 150 IPARD = 1 151 IF (MR.NE.0) IPARD = 0 152 INITD = 0 153 CALL ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, PHIDR, PHIDI, 154 * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, CWRKR(1,3), 155 * CWRKI(1,3)) 156 IF (KODE.EQ.1) GO TO 80 157 STR = ZRR + ZET2DR 158 STI = ZRI + ZET2DI 159 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 160 STR = STR*RAST*RAST 161 STI = -STI*RAST*RAST 162 S1R = ZET1DR - STR 163 S1I = ZET1DI - STI 164 GO TO 90 165 80 CONTINUE 166 S1R = ZET1DR - ZET2DR 167 S1I = ZET1DI - ZET2DI 168 90 CONTINUE 169 RS1 = S1R 170 IF (DABS(RS1).GT.ELIM) GO TO 95 171 IF (DABS(RS1).LT.ALIM) GO TO 100 172 C---------------------------------------------------------------------------- 173 C REFINE ESTIMATE AND TEST 174 C------------------------------------------------------------------------- 175 APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0))) 176 RS1 = RS1+DLOG(APHI) 177 IF (DABS(RS1).LT.ELIM) GO TO 100 178 95 CONTINUE 179 IF (DABS(RS1).GT.0.0D0) GO TO 300 180 C----------------------------------------------------------------------- 181 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 182 C----------------------------------------------------------------------- 183 IF (ZR.LT.0.0D0) GO TO 300 184 NZ = N 185 DO 96 I=1,N 186 YR(I) = ZEROR 187 YI(I) = ZEROI 188 96 CONTINUE 189 RETURN 190 C--------------------------------------------------------------------------- 191 C FORWARD RECUR FOR REMAINDER OF THE SEQUENCE 192 C---------------------------------------------------------------------------- 193 100 CONTINUE 194 S1R = CYR(1) 195 S1I = CYI(1) 196 S2R = CYR(2) 197 S2I = CYI(2) 198 C1R = CSRR(KFLAG) 199 ASCLE = BRY(KFLAG) 200 DO 120 I=IB,N 201 C2R = S2R 202 C2I = S2I 203 S2R = CKR*C2R - CKI*C2I + S1R 204 S2I = CKR*C2I + CKI*C2R + S1I 205 S1R = C2R 206 S1I = C2I 207 CKR = CKR + RZR 208 CKI = CKI + RZI 209 C2R = S2R*C1R 210 C2I = S2I*C1R 211 YR(I) = C2R 212 YI(I) = C2I 213 IF (KFLAG.GE.3) GO TO 120 214 STR = DABS(C2R) 215 STI = DABS(C2I) 216 C2M = DMAX1(STR,STI) 217 IF (C2M.LE.ASCLE) GO TO 120 218 KFLAG = KFLAG + 1 219 ASCLE = BRY(KFLAG) 220 S1R = S1R*C1R 221 S1I = S1I*C1R 222 S2R = C2R 223 S2I = C2I 224 S1R = S1R*CSSR(KFLAG) 225 S1I = S1I*CSSR(KFLAG) 226 S2R = S2R*CSSR(KFLAG) 227 S2I = S2I*CSSR(KFLAG) 228 C1R = CSRR(KFLAG) 229 120 CONTINUE 230 160 CONTINUE 231 IF (MR.EQ.0) RETURN 232 C----------------------------------------------------------------------- 233 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0 234 C----------------------------------------------------------------------- 235 NZ = 0 236 FMR = DBLE(FLOAT(MR)) 237 SGN = -DSIGN(PI,FMR) 238 C----------------------------------------------------------------------- 239 C CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP. 240 C----------------------------------------------------------------------- 241 CSGNI = SGN 242 INU = INT(SNGL(FNU)) 243 FNF = FNU - DBLE(FLOAT(INU)) 244 IFN = INU + N - 1 245 ANG = FNF*SGN 246 CSPNR = DCOS(ANG) 247 CSPNI = DSIN(ANG) 248 IF (MOD(IFN,2).EQ.0) GO TO 170 249 CSPNR = -CSPNR 250 CSPNI = -CSPNI 251 170 CONTINUE 252 ASC = BRY(1) 253 IUF = 0 254 KK = N 255 KDFLG = 1 256 IB = IB - 1 257 IC = IB - 1 258 DO 270 K=1,N 259 FN = FNU + DBLE(FLOAT(KK-1)) 260 C----------------------------------------------------------------------- 261 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K 262 C FUNCTION ABOVE 263 C----------------------------------------------------------------------- 264 M=3 265 IF (N.GT.2) GO TO 175 266 172 CONTINUE 267 INITD = INIT(J) 268 PHIDR = PHIR(J) 269 PHIDI = PHII(J) 270 ZET1DR = ZETA1R(J) 271 ZET1DI = ZETA1I(J) 272 ZET2DR = ZETA2R(J) 273 ZET2DI = ZETA2I(J) 274 SUMDR = SUMR(J) 275 SUMDI = SUMI(J) 276 M = J 277 J = 3 - J 278 GO TO 180 279 175 CONTINUE 280 IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180 281 IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172 282 INITD = 0 283 180 CONTINUE 284 CALL ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, PHIDR, PHIDI, 285 * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, 286 * CWRKR(1,M), CWRKI(1,M)) 287 IF (KODE.EQ.1) GO TO 200 288 STR = ZRR + ZET2DR 289 STI = ZRI + ZET2DI 290 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 291 STR = STR*RAST*RAST 292 STI = -STI*RAST*RAST 293 S1R = -ZET1DR + STR 294 S1I = -ZET1DI + STI 295 GO TO 210 296 200 CONTINUE 297 S1R = -ZET1DR + ZET2DR 298 S1I = -ZET1DI + ZET2DI 299 210 CONTINUE 300 C----------------------------------------------------------------------- 301 C TEST FOR UNDERFLOW AND OVERFLOW 302 C----------------------------------------------------------------------- 303 RS1 = S1R 304 IF (DABS(RS1).GT.ELIM) GO TO 260 305 IF (KDFLG.EQ.1) IFLAG = 2 306 IF (DABS(RS1).LT.ALIM) GO TO 220 307 C----------------------------------------------------------------------- 308 C REFINE TEST AND SCALE 309 C----------------------------------------------------------------------- 310 APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0))) 311 RS1 = RS1 + DLOG(APHI) 312 IF (DABS(RS1).GT.ELIM) GO TO 260 313 IF (KDFLG.EQ.1) IFLAG = 1 314 IF (RS1.LT.0.0D0) GO TO 220 315 IF (KDFLG.EQ.1) IFLAG = 3 316 220 CONTINUE 317 STR = PHIDR*SUMDR - PHIDI*SUMDI 318 STI = PHIDR*SUMDI + PHIDI*SUMDR 319 S2R = -CSGNI*STI 320 S2I = CSGNI*STR 321 STR = DEXP(S1R)*CSSR(IFLAG) 322 S1R = STR*DCOS(S1I) 323 S1I = STR*DSIN(S1I) 324 STR = S2R*S1R - S2I*S1I 325 S2I = S2R*S1I + S2I*S1R 326 S2R = STR 327 IF (IFLAG.NE.1) GO TO 230 328 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 329 IF (NW.EQ.0) GO TO 230 330 S2R = ZEROR 331 S2I = ZEROI 332 230 CONTINUE 333 CYR(KDFLG) = S2R 334 CYI(KDFLG) = S2I 335 C2R = S2R 336 C2I = S2I 337 S2R = S2R*CSRR(IFLAG) 338 S2I = S2I*CSRR(IFLAG) 339 C----------------------------------------------------------------------- 340 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N 341 C----------------------------------------------------------------------- 342 S1R = YR(KK) 343 S1I = YI(KK) 344 IF (KODE.EQ.1) GO TO 250 345 CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF) 346 NZ = NZ + NW 347 250 CONTINUE 348 YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R 349 YI(KK) = CSPNR*S1I + CSPNI*S1R + S2I 350 KK = KK - 1 351 CSPNR = -CSPNR 352 CSPNI = -CSPNI 353 IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255 354 KDFLG = 1 355 GO TO 270 356 255 CONTINUE 357 IF (KDFLG.EQ.2) GO TO 275 358 KDFLG = 2 359 GO TO 270 360 260 CONTINUE 361 IF (RS1.GT.0.0D0) GO TO 300 362 S2R = ZEROR 363 S2I = ZEROI 364 GO TO 230 365 270 CONTINUE 366 K = N 367 275 CONTINUE 368 IL = N - K 369 IF (IL.EQ.0) RETURN 370 C----------------------------------------------------------------------- 371 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE 372 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP 373 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES. 374 C----------------------------------------------------------------------- 375 S1R = CYR(1) 376 S1I = CYI(1) 377 S2R = CYR(2) 378 S2I = CYI(2) 379 CSR = CSRR(IFLAG) 380 ASCLE = BRY(IFLAG) 381 FN = DBLE(FLOAT(INU+IL)) 382 DO 290 I=1,IL 383 C2R = S2R 384 C2I = S2I 385 S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I) 386 S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R) 387 S1R = C2R 388 S1I = C2I 389 FN = FN - 1.0D0 390 C2R = S2R*CSR 391 C2I = S2I*CSR 392 CKR = C2R 393 CKI = C2I 394 C1R = YR(KK) 395 C1I = YI(KK) 396 IF (KODE.EQ.1) GO TO 280 397 CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF) 398 NZ = NZ + NW 399 280 CONTINUE 400 YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R 401 YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I 402 KK = KK - 1 403 CSPNR = -CSPNR 404 CSPNI = -CSPNI 405 IF (IFLAG.GE.3) GO TO 290 406 C2R = DABS(CKR) 407 C2I = DABS(CKI) 408 C2M = DMAX1(C2R,C2I) 409 IF (C2M.LE.ASCLE) GO TO 290 410 IFLAG = IFLAG + 1 411 ASCLE = BRY(IFLAG) 412 S1R = S1R*CSR 413 S1I = S1I*CSR 414 S2R = CKR 415 S2I = CKI 416 S1R = S1R*CSSR(IFLAG) 417 S1I = S1I*CSSR(IFLAG) 418 S2R = S2R*CSSR(IFLAG) 419 S2I = S2I*CSSR(IFLAG) 420 CSR = CSRR(IFLAG) 421 290 CONTINUE 422 RETURN 423 300 CONTINUE 424 NZ = -1 425 RETURN 426 END