src/numerical/slatec/fortran/zunk2.f

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