github.com/gopherd/gonum@v0.0.4/mathext/internal/amos/amoslib/zunk2.f (about) 1 SUBROUTINE ZUNK2(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM, 2 * ALIM) 3 C***BEGIN PROLOGUE ZUNK2 4 C***REFER TO ZBESK 5 C 6 C ZUNK2 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 EXPANSIONS FOR H(KIND,FNU,ZN) AND J(FNU,ZN) 9 C WHERE ZN IS IN THE RIGHT HALF PLANE, KIND=(3-MR)/2, MR=+1 OR 10 C -1. HERE ZN=ZR*I OR -ZR*I WHERE ZR=Z IF Z IS IN THE RIGHT 11 C HALF PLANE OR ZR=-Z IF Z IS IN THE LEFT HALF PLANE. MR INDIC- 12 C ATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION. 13 C NZ=-1 MEANS AN OVERFLOW WILL OCCUR 14 C 15 C***ROUTINES CALLED ZAIRY,ZKSCL,ZS1S2,ZUCHK,ZUNHJ,D1MACH,ZABS 16 C***END PROLOGUE ZUNK2 17 C COMPLEX AI,ARG,ARGD,ASUM,ASUMD,BSUM,BSUMD,CFN,CI,CIP,CK,CONE,CRSC, 18 C *CR1,CR2,CS,CSCL,CSGN,CSPN,CSR,CSS,CY,CZERO,C1,C2,DAI,PHI,PHID,RZ, 19 C *S1,S2,Y,Z,ZB,ZETA1,ZETA1D,ZETA2,ZETA2D,ZN,ZR 20 DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGDI, 21 * ARGDR, ARGI, ARGR, ASC, ASCLE, ASUMDI, ASUMDR, ASUMI, ASUMR, 22 * BRY, BSUMDI, BSUMDR, BSUMI, BSUMR, CAR, CIPI, CIPR, CKI, CKR, 23 * CONER, CRSC, CR1I, CR1R, CR2I, CR2R, CSCL, CSGNI, CSI, 24 * CSPNI, CSPNR, CSR, CSRR, CSSR, CYI, CYR, C1I, C1R, C2I, C2M, 25 * C2R, DAII, DAIR, ELIM, FMR, FN, FNF, FNU, HPI, PHIDI, PHIDR, 26 * PHII, PHIR, PI, PTI, PTR, RAST, RAZR, RS1, RZI, RZR, SAR, SGN, 27 * STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, YY, ZBI, ZBR, ZEROI, 28 * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZET1DI, ZET1DR, ZET2DI, 29 * ZET2DR, ZI, ZNI, ZNR, ZR, ZRI, ZRR, D1MACH, ZABS 30 INTEGER I, IB, IFLAG, IFN, IL, IN, INU, IUF, K, KDFLG, KFLAG, KK, 31 * KODE, MR, N, NAI, NDAI, NW, NZ, IDUM, J, IPARD, IC 32 DIMENSION BRY(3), YR(N), YI(N), ASUMR(2), ASUMI(2), BSUMR(2), 33 * BSUMI(2), PHIR(2), PHII(2), ARGR(2), ARGI(2), ZETA1R(2), 34 * ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2), CIPR(4), 35 * CIPI(4), CSSR(3), CSRR(3) 36 DATA ZEROR,ZEROI,CONER,CR1R,CR1I,CR2R,CR2I / 37 1 0.0D0, 0.0D0, 1.0D0, 38 1 1.0D0,1.73205080756887729D0 , -0.5D0,-8.66025403784438647D-01 / 39 DATA HPI, PI, AIC / 40 1 1.57079632679489662D+00, 3.14159265358979324D+00, 41 1 1.26551212348464539D+00/ 42 DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4), 43 * CIPI(4) / 44 1 1.0D0,0.0D0 , 0.0D0,-1.0D0 , -1.0D0,0.0D0 , 0.0D0,1.0D0 / 45 C 46 KDFLG = 1 47 NZ = 0 48 C----------------------------------------------------------------------- 49 C EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN 50 C THE UNDERFLOW LIMIT 51 C----------------------------------------------------------------------- 52 CSCL = 1.0D0/TOL 53 CRSC = TOL 54 CSSR(1) = CSCL 55 CSSR(2) = CONER 56 CSSR(3) = CRSC 57 CSRR(1) = CRSC 58 CSRR(2) = CONER 59 CSRR(3) = CSCL 60 BRY(1) = 1.0D+3*D1MACH(1)/TOL 61 BRY(2) = 1.0D0/BRY(1) 62 BRY(3) = D1MACH(2) 63 ZRR = ZR 64 ZRI = ZI 65 IF (ZR.GE.0.0D0) GO TO 10 66 ZRR = -ZR 67 ZRI = -ZI 68 10 CONTINUE 69 YY = ZRI 70 ZNR = ZRI 71 ZNI = -ZRR 72 ZBR = ZRR 73 ZBI = ZRI 74 INU = INT(SNGL(FNU)) 75 FNF = FNU - DBLE(FLOAT(INU)) 76 ANG = -HPI*FNF 77 CAR = DCOS(ANG) 78 SAR = DSIN(ANG) 79 C2R = HPI*SAR 80 C2I = -HPI*CAR 81 KK = MOD(INU,4) + 1 82 STR = C2R*CIPR(KK) - C2I*CIPI(KK) 83 STI = C2R*CIPI(KK) + C2I*CIPR(KK) 84 CSR = CR1R*STR - CR1I*STI 85 CSI = CR1R*STI + CR1I*STR 86 IF (YY.GT.0.0D0) GO TO 20 87 ZNR = -ZNR 88 ZBI = -ZBI 89 20 CONTINUE 90 C----------------------------------------------------------------------- 91 C K(FNU,Z) IS COMPUTED FROM H(2,FNU,-I*Z) WHERE Z IS IN THE FIRST 92 C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY 93 C CONJUGATION SINCE THE K FUNCTION IS REAL ON THE POSITIVE REAL AXIS 94 C----------------------------------------------------------------------- 95 J = 2 96 DO 80 I=1,N 97 C----------------------------------------------------------------------- 98 C J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J 99 C----------------------------------------------------------------------- 100 J = 3 - J 101 FN = FNU + DBLE(FLOAT(I-1)) 102 CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR(J), PHII(J), ARGR(J), 103 * ARGI(J), ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), ASUMR(J), 104 * ASUMI(J), BSUMR(J), BSUMI(J)) 105 IF (KODE.EQ.1) GO TO 30 106 STR = ZBR + ZETA2R(J) 107 STI = ZBI + ZETA2I(J) 108 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 109 STR = STR*RAST*RAST 110 STI = -STI*RAST*RAST 111 S1R = ZETA1R(J) - STR 112 S1I = ZETA1I(J) - STI 113 GO TO 40 114 30 CONTINUE 115 S1R = ZETA1R(J) - ZETA2R(J) 116 S1I = ZETA1I(J) - ZETA2I(J) 117 40 CONTINUE 118 C----------------------------------------------------------------------- 119 C TEST FOR UNDERFLOW AND OVERFLOW 120 C----------------------------------------------------------------------- 121 RS1 = S1R 122 IF (DABS(RS1).GT.ELIM) GO TO 70 123 IF (KDFLG.EQ.1) KFLAG = 2 124 IF (DABS(RS1).LT.ALIM) GO TO 50 125 C----------------------------------------------------------------------- 126 C REFINE TEST AND SCALE 127 C----------------------------------------------------------------------- 128 APHI = ZABS(CMPLX(PHIR(J),PHII(J),kind=KIND(1.0D0))) 129 AARG = ZABS(CMPLX(ARGR(J),ARGI(J),kind=KIND(1.0D0))) 130 RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC 131 IF (DABS(RS1).GT.ELIM) GO TO 70 132 IF (KDFLG.EQ.1) KFLAG = 1 133 IF (RS1.LT.0.0D0) GO TO 50 134 IF (KDFLG.EQ.1) KFLAG = 3 135 50 CONTINUE 136 C----------------------------------------------------------------------- 137 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR 138 C EXPONENT EXTREMES 139 C----------------------------------------------------------------------- 140 C2R = ARGR(J)*CR2R - ARGI(J)*CR2I 141 C2I = ARGR(J)*CR2I + ARGI(J)*CR2R 142 CALL ZAIRY(C2R, C2I, 0, 2, AIR, AII, NAI, IDUM) 143 CALL ZAIRY(C2R, C2I, 1, 2, DAIR, DAII, NDAI, IDUM) 144 STR = DAIR*BSUMR(J) - DAII*BSUMI(J) 145 STI = DAIR*BSUMI(J) + DAII*BSUMR(J) 146 PTR = STR*CR2R - STI*CR2I 147 PTI = STR*CR2I + STI*CR2R 148 STR = PTR + (AIR*ASUMR(J)-AII*ASUMI(J)) 149 STI = PTI + (AIR*ASUMI(J)+AII*ASUMR(J)) 150 PTR = STR*PHIR(J) - STI*PHII(J) 151 PTI = STR*PHII(J) + STI*PHIR(J) 152 S2R = PTR*CSR - PTI*CSI 153 S2I = PTR*CSI + PTI*CSR 154 STR = DEXP(S1R)*CSSR(KFLAG) 155 S1R = STR*DCOS(S1I) 156 S1I = STR*DSIN(S1I) 157 STR = S2R*S1R - S2I*S1I 158 S2I = S1R*S2I + S2R*S1I 159 S2R = STR 160 IF (KFLAG.NE.1) GO TO 60 161 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 162 IF (NW.NE.0) GO TO 70 163 60 CONTINUE 164 IF (YY.LE.0.0D0) S2I = -S2I 165 CYR(KDFLG) = S2R 166 CYI(KDFLG) = S2I 167 YR(I) = S2R*CSRR(KFLAG) 168 YI(I) = S2I*CSRR(KFLAG) 169 STR = CSI 170 CSI = -CSR 171 CSR = STR 172 IF (KDFLG.EQ.2) GO TO 85 173 KDFLG = 2 174 GO TO 80 175 70 CONTINUE 176 IF (RS1.GT.0.0D0) GO TO 320 177 C----------------------------------------------------------------------- 178 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 179 C----------------------------------------------------------------------- 180 IF (ZR.LT.0.0D0) GO TO 320 181 KDFLG = 1 182 YR(I)=ZEROR 183 YI(I)=ZEROI 184 NZ=NZ+1 185 STR = CSI 186 CSI =-CSR 187 CSR = STR 188 IF (I.EQ.1) GO TO 80 189 IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 80 190 YR(I-1)=ZEROR 191 YI(I-1)=ZEROI 192 NZ=NZ+1 193 80 CONTINUE 194 I = N 195 85 CONTINUE 196 RAZR = 1.0D0/ZABS(CMPLX(ZRR,ZRI,kind=KIND(1.0D0))) 197 STR = ZRR*RAZR 198 STI = -ZRI*RAZR 199 RZR = (STR+STR)*RAZR 200 RZI = (STI+STI)*RAZR 201 CKR = FN*RZR 202 CKI = FN*RZI 203 IB = I + 1 204 IF (N.LT.IB) GO TO 180 205 C----------------------------------------------------------------------- 206 C TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO 207 C ON UNDERFLOW. 208 C----------------------------------------------------------------------- 209 FN = FNU + DBLE(FLOAT(N-1)) 210 IPARD = 1 211 IF (MR.NE.0) IPARD = 0 212 CALL ZUNHJ(ZNR, ZNI, FN, IPARD, TOL, PHIDR, PHIDI, ARGDR, ARGDI, 213 * ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR, ASUMDI, BSUMDR, BSUMDI) 214 IF (KODE.EQ.1) GO TO 90 215 STR = ZBR + ZET2DR 216 STI = ZBI + ZET2DI 217 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 218 STR = STR*RAST*RAST 219 STI = -STI*RAST*RAST 220 S1R = ZET1DR - STR 221 S1I = ZET1DI - STI 222 GO TO 100 223 90 CONTINUE 224 S1R = ZET1DR - ZET2DR 225 S1I = ZET1DI - ZET2DI 226 100 CONTINUE 227 RS1 = S1R 228 IF (DABS(RS1).GT.ELIM) GO TO 105 229 IF (DABS(RS1).LT.ALIM) GO TO 120 230 C---------------------------------------------------------------------------- 231 C REFINE ESTIMATE AND TEST 232 C------------------------------------------------------------------------- 233 APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0))) 234 RS1 = RS1+DLOG(APHI) 235 IF (DABS(RS1).LT.ELIM) GO TO 120 236 105 CONTINUE 237 IF (RS1.GT.0.0D0) GO TO 320 238 C----------------------------------------------------------------------- 239 C FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW 240 C----------------------------------------------------------------------- 241 IF (ZR.LT.0.0D0) GO TO 320 242 NZ = N 243 DO 106 I=1,N 244 YR(I) = ZEROR 245 YI(I) = ZEROI 246 106 CONTINUE 247 RETURN 248 120 CONTINUE 249 S1R = CYR(1) 250 S1I = CYI(1) 251 S2R = CYR(2) 252 S2I = CYI(2) 253 C1R = CSRR(KFLAG) 254 ASCLE = BRY(KFLAG) 255 DO 130 I=IB,N 256 C2R = S2R 257 C2I = S2I 258 S2R = CKR*C2R - CKI*C2I + S1R 259 S2I = CKR*C2I + CKI*C2R + S1I 260 S1R = C2R 261 S1I = C2I 262 CKR = CKR + RZR 263 CKI = CKI + RZI 264 C2R = S2R*C1R 265 C2I = S2I*C1R 266 YR(I) = C2R 267 YI(I) = C2I 268 IF (KFLAG.GE.3) GO TO 130 269 STR = DABS(C2R) 270 STI = DABS(C2I) 271 C2M = DMAX1(STR,STI) 272 IF (C2M.LE.ASCLE) GO TO 130 273 KFLAG = KFLAG + 1 274 ASCLE = BRY(KFLAG) 275 S1R = S1R*C1R 276 S1I = S1I*C1R 277 S2R = C2R 278 S2I = C2I 279 S1R = S1R*CSSR(KFLAG) 280 S1I = S1I*CSSR(KFLAG) 281 S2R = S2R*CSSR(KFLAG) 282 S2I = S2I*CSSR(KFLAG) 283 C1R = CSRR(KFLAG) 284 130 CONTINUE 285 180 CONTINUE 286 IF (MR.EQ.0) RETURN 287 C----------------------------------------------------------------------- 288 C ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0 289 C----------------------------------------------------------------------- 290 NZ = 0 291 FMR = DBLE(FLOAT(MR)) 292 SGN = -DSIGN(PI,FMR) 293 C----------------------------------------------------------------------- 294 C CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP. 295 C----------------------------------------------------------------------- 296 CSGNI = SGN 297 IF (YY.LE.0.0D0) CSGNI = -CSGNI 298 IFN = INU + N - 1 299 ANG = FNF*SGN 300 CSPNR = DCOS(ANG) 301 CSPNI = DSIN(ANG) 302 IF (MOD(IFN,2).EQ.0) GO TO 190 303 CSPNR = -CSPNR 304 CSPNI = -CSPNI 305 190 CONTINUE 306 C----------------------------------------------------------------------- 307 C CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS 308 C COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST 309 C QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY 310 C CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS 311 C----------------------------------------------------------------------- 312 CSR = SAR*CSGNI 313 CSI = CAR*CSGNI 314 IN = MOD(IFN,4) + 1 315 C2R = CIPR(IN) 316 C2I = CIPI(IN) 317 STR = CSR*C2R + CSI*C2I 318 CSI = -CSR*C2I + CSI*C2R 319 CSR = STR 320 ASC = BRY(1) 321 IUF = 0 322 KK = N 323 KDFLG = 1 324 IB = IB - 1 325 IC = IB - 1 326 DO 290 K=1,N 327 FN = FNU + DBLE(FLOAT(KK-1)) 328 C----------------------------------------------------------------------- 329 C LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K 330 C FUNCTION ABOVE 331 C----------------------------------------------------------------------- 332 IF (N.GT.2) GO TO 175 333 172 CONTINUE 334 PHIDR = PHIR(J) 335 PHIDI = PHII(J) 336 ARGDR = ARGR(J) 337 ARGDI = ARGI(J) 338 ZET1DR = ZETA1R(J) 339 ZET1DI = ZETA1I(J) 340 ZET2DR = ZETA2R(J) 341 ZET2DI = ZETA2I(J) 342 ASUMDR = ASUMR(J) 343 ASUMDI = ASUMI(J) 344 BSUMDR = BSUMR(J) 345 BSUMDI = BSUMI(J) 346 J = 3 - J 347 GO TO 210 348 175 CONTINUE 349 IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 210 350 IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172 351 CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIDR, PHIDI, ARGDR, 352 * ARGDI, ZET1DR, ZET1DI, ZET2DR, ZET2DI, ASUMDR, 353 * ASUMDI, BSUMDR, BSUMDI) 354 210 CONTINUE 355 IF (KODE.EQ.1) GO TO 220 356 STR = ZBR + ZET2DR 357 STI = ZBI + ZET2DI 358 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 359 STR = STR*RAST*RAST 360 STI = -STI*RAST*RAST 361 S1R = -ZET1DR + STR 362 S1I = -ZET1DI + STI 363 GO TO 230 364 220 CONTINUE 365 S1R = -ZET1DR + ZET2DR 366 S1I = -ZET1DI + ZET2DI 367 230 CONTINUE 368 C----------------------------------------------------------------------- 369 C TEST FOR UNDERFLOW AND OVERFLOW 370 C----------------------------------------------------------------------- 371 RS1 = S1R 372 IF (DABS(RS1).GT.ELIM) GO TO 280 373 IF (KDFLG.EQ.1) IFLAG = 2 374 IF (DABS(RS1).LT.ALIM) GO TO 240 375 C----------------------------------------------------------------------- 376 C REFINE TEST AND SCALE 377 C----------------------------------------------------------------------- 378 APHI = ZABS(CMPLX(PHIDR,PHIDI,kind=KIND(1.0D0))) 379 AARG = ZABS(CMPLX(ARGDR,ARGDI,kind=KIND(1.0D0))) 380 RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC 381 IF (DABS(RS1).GT.ELIM) GO TO 280 382 IF (KDFLG.EQ.1) IFLAG = 1 383 IF (RS1.LT.0.0D0) GO TO 240 384 IF (KDFLG.EQ.1) IFLAG = 3 385 240 CONTINUE 386 CALL ZAIRY(ARGDR, ARGDI, 0, 2, AIR, AII, NAI, IDUM) 387 CALL ZAIRY(ARGDR, ARGDI, 1, 2, DAIR, DAII, NDAI, IDUM) 388 STR = DAIR*BSUMDR - DAII*BSUMDI 389 STI = DAIR*BSUMDI + DAII*BSUMDR 390 STR = STR + (AIR*ASUMDR-AII*ASUMDI) 391 STI = STI + (AIR*ASUMDI+AII*ASUMDR) 392 PTR = STR*PHIDR - STI*PHIDI 393 PTI = STR*PHIDI + STI*PHIDR 394 S2R = PTR*CSR - PTI*CSI 395 S2I = PTR*CSI + PTI*CSR 396 STR = DEXP(S1R)*CSSR(IFLAG) 397 S1R = STR*DCOS(S1I) 398 S1I = STR*DSIN(S1I) 399 STR = S2R*S1R - S2I*S1I 400 S2I = S2R*S1I + S2I*S1R 401 S2R = STR 402 IF (IFLAG.NE.1) GO TO 250 403 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 404 IF (NW.EQ.0) GO TO 250 405 S2R = ZEROR 406 S2I = ZEROI 407 250 CONTINUE 408 IF (YY.LE.0.0D0) S2I = -S2I 409 CYR(KDFLG) = S2R 410 CYI(KDFLG) = S2I 411 C2R = S2R 412 C2I = S2I 413 S2R = S2R*CSRR(IFLAG) 414 S2I = S2I*CSRR(IFLAG) 415 C----------------------------------------------------------------------- 416 C ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N 417 C----------------------------------------------------------------------- 418 S1R = YR(KK) 419 S1I = YI(KK) 420 IF (KODE.EQ.1) GO TO 270 421 CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF) 422 NZ = NZ + NW 423 270 CONTINUE 424 YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R 425 YI(KK) = S1R*CSPNI + S1I*CSPNR + S2I 426 KK = KK - 1 427 CSPNR = -CSPNR 428 CSPNI = -CSPNI 429 STR = CSI 430 CSI = -CSR 431 CSR = STR 432 IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255 433 KDFLG = 1 434 GO TO 290 435 255 CONTINUE 436 IF (KDFLG.EQ.2) GO TO 295 437 KDFLG = 2 438 GO TO 290 439 280 CONTINUE 440 IF (RS1.GT.0.0D0) GO TO 320 441 S2R = ZEROR 442 S2I = ZEROI 443 GO TO 250 444 290 CONTINUE 445 K = N 446 295 CONTINUE 447 IL = N - K 448 IF (IL.EQ.0) RETURN 449 C----------------------------------------------------------------------- 450 C RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE 451 C K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP 452 C INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES. 453 C----------------------------------------------------------------------- 454 S1R = CYR(1) 455 S1I = CYI(1) 456 S2R = CYR(2) 457 S2I = CYI(2) 458 CSR = CSRR(IFLAG) 459 ASCLE = BRY(IFLAG) 460 FN = DBLE(FLOAT(INU+IL)) 461 DO 310 I=1,IL 462 C2R = S2R 463 C2I = S2I 464 S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I) 465 S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R) 466 S1R = C2R 467 S1I = C2I 468 FN = FN - 1.0D0 469 C2R = S2R*CSR 470 C2I = S2I*CSR 471 CKR = C2R 472 CKI = C2I 473 C1R = YR(KK) 474 C1I = YI(KK) 475 IF (KODE.EQ.1) GO TO 300 476 CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF) 477 NZ = NZ + NW 478 300 CONTINUE 479 YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R 480 YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I 481 KK = KK - 1 482 CSPNR = -CSPNR 483 CSPNI = -CSPNI 484 IF (IFLAG.GE.3) GO TO 310 485 C2R = DABS(CKR) 486 C2I = DABS(CKI) 487 C2M = DMAX1(C2R,C2I) 488 IF (C2M.LE.ASCLE) GO TO 310 489 IFLAG = IFLAG + 1 490 ASCLE = BRY(IFLAG) 491 S1R = S1R*CSR 492 S1I = S1I*CSR 493 S2R = CKR 494 S2I = CKI 495 S1R = S1R*CSSR(IFLAG) 496 S1I = S1I*CSSR(IFLAG) 497 S2R = S2R*CSSR(IFLAG) 498 S2I = S2I*CSSR(IFLAG) 499 CSR = CSRR(IFLAG) 500 310 CONTINUE 501 RETURN 502 320 CONTINUE 503 NZ = -1 504 RETURN 505 END