github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zuni2.f (about) 1 SUBROUTINE ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL, 2 * TOL, ELIM, ALIM) 3 C***BEGIN PROLOGUE ZUNI2 4 C***REFER TO ZBESI,ZBESK 5 C 6 C ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF 7 C UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I 8 C OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO. 9 C 10 C FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC 11 C EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET. 12 C NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER 13 C FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL. 14 C Y(I)=CZERO FOR I=NLAST+1,N 15 C 16 C***ROUTINES CALLED ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,ZABS 17 C***END PROLOGUE ZUNI2 18 C COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS, 19 C *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN 20 DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI, 21 * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR, 22 * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII, 23 * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI, 24 * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI, 25 * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR, 26 * CYI, D1MACH, ZABS, CAR, SAR 27 INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST, 28 * NN, NUF, NW, NZ, IDUM 29 DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3), 30 * CSRR(3), CYR(2), CYI(2) 31 DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 / 32 DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4), 33 * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/ 34 DATA HPI, AIC / 35 1 1.57079632679489662D+00, 1.265512123484645396D+00/ 36 C 37 NZ = 0 38 ND = N 39 NLAST = 0 40 C----------------------------------------------------------------------- 41 C COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG- 42 C NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE, 43 C EXP(ALIM)=EXP(ELIM)*TOL 44 C----------------------------------------------------------------------- 45 CSCL = 1.0D0/TOL 46 CRSC = TOL 47 CSSR(1) = CSCL 48 CSSR(2) = CONER 49 CSSR(3) = CRSC 50 CSRR(1) = CRSC 51 CSRR(2) = CONER 52 CSRR(3) = CSCL 53 BRY(1) = 1.0D+3*D1MACH(1)/TOL 54 C----------------------------------------------------------------------- 55 C ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI 56 C----------------------------------------------------------------------- 57 ZNR = ZI 58 ZNI = -ZR 59 ZBR = ZR 60 ZBI = ZI 61 CIDI = -CONER 62 INU = INT(SNGL(FNU)) 63 ANG = HPI*(FNU-DBLE(FLOAT(INU))) 64 C2R = DCOS(ANG) 65 C2I = DSIN(ANG) 66 CAR = C2R 67 SAR = C2I 68 IN = INU + N - 1 69 IN = MOD(IN,4) + 1 70 STR = C2R*CIPR(IN) - C2I*CIPI(IN) 71 C2I = C2R*CIPI(IN) + C2I*CIPR(IN) 72 C2R = STR 73 IF (ZI.GT.0.0D0) GO TO 10 74 ZNR = -ZNR 75 ZBI = -ZBI 76 CIDI = -CIDI 77 C2I = -C2I 78 10 CONTINUE 79 C----------------------------------------------------------------------- 80 C CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER 81 C----------------------------------------------------------------------- 82 FN = DMAX1(FNU,1.0D0) 83 CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R, 84 * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) 85 IF (KODE.EQ.1) GO TO 20 86 STR = ZBR + ZETA2R 87 STI = ZBI + ZETA2I 88 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 89 STR = STR*RAST*RAST 90 STI = -STI*RAST*RAST 91 S1R = -ZETA1R + STR 92 S1I = -ZETA1I + STI 93 GO TO 30 94 20 CONTINUE 95 S1R = -ZETA1R + ZETA2R 96 S1I = -ZETA1I + ZETA2I 97 30 CONTINUE 98 RS1 = S1R 99 IF (DABS(RS1).GT.ELIM) GO TO 150 100 40 CONTINUE 101 NN = MIN0(2,ND) 102 DO 90 I=1,NN 103 FN = FNU + DBLE(FLOAT(ND-I)) 104 CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI, 105 * ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI) 106 IF (KODE.EQ.1) GO TO 50 107 STR = ZBR + ZETA2R 108 STI = ZBI + ZETA2I 109 RAST = FN/ZABS(CMPLX(STR,STI,kind=KIND(1.0D0))) 110 STR = STR*RAST*RAST 111 STI = -STI*RAST*RAST 112 S1R = -ZETA1R + STR 113 S1I = -ZETA1I + STI + DABS(ZI) 114 GO TO 60 115 50 CONTINUE 116 S1R = -ZETA1R + ZETA2R 117 S1I = -ZETA1I + ZETA2I 118 60 CONTINUE 119 C----------------------------------------------------------------------- 120 C TEST FOR UNDERFLOW AND OVERFLOW 121 C----------------------------------------------------------------------- 122 RS1 = S1R 123 IF (DABS(RS1).GT.ELIM) GO TO 120 124 IF (I.EQ.1) IFLAG = 2 125 IF (DABS(RS1).LT.ALIM) GO TO 70 126 C----------------------------------------------------------------------- 127 C REFINE TEST AND SCALE 128 C----------------------------------------------------------------------- 129 C----------------------------------------------------------------------- 130 APHI = ZABS(CMPLX(PHIR,PHII,kind=KIND(1.0D0))) 131 AARG = ZABS(CMPLX(ARGR,ARGI,kind=KIND(1.0D0))) 132 RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC 133 IF (DABS(RS1).GT.ELIM) GO TO 120 134 IF (I.EQ.1) IFLAG = 1 135 IF (RS1.LT.0.0D0) GO TO 70 136 IF (I.EQ.1) IFLAG = 3 137 70 CONTINUE 138 C----------------------------------------------------------------------- 139 C SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR 140 C EXPONENT EXTREMES 141 C----------------------------------------------------------------------- 142 CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM) 143 CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM) 144 STR = DAIR*BSUMR - DAII*BSUMI 145 STI = DAIR*BSUMI + DAII*BSUMR 146 STR = STR + (AIR*ASUMR-AII*ASUMI) 147 STI = STI + (AIR*ASUMI+AII*ASUMR) 148 S2R = PHIR*STR - PHII*STI 149 S2I = PHIR*STI + PHII*STR 150 STR = DEXP(S1R)*CSSR(IFLAG) 151 S1R = STR*DCOS(S1I) 152 S1I = STR*DSIN(S1I) 153 STR = S2R*S1R - S2I*S1I 154 S2I = S2R*S1I + S2I*S1R 155 S2R = STR 156 IF (IFLAG.NE.1) GO TO 80 157 CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL) 158 IF (NW.NE.0) GO TO 120 159 80 CONTINUE 160 IF (ZI.LE.0.0D0) S2I = -S2I 161 STR = S2R*C2R - S2I*C2I 162 S2I = S2R*C2I + S2I*C2R 163 S2R = STR 164 CYR(I) = S2R 165 CYI(I) = S2I 166 J = ND - I + 1 167 YR(J) = S2R*CSRR(IFLAG) 168 YI(J) = S2I*CSRR(IFLAG) 169 STR = -C2I*CIDI 170 C2I = C2R*CIDI 171 C2R = STR 172 90 CONTINUE 173 IF (ND.LE.2) GO TO 110 174 RAZ = 1.0D0/ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 175 STR = ZR*RAZ 176 STI = -ZI*RAZ 177 RZR = (STR+STR)*RAZ 178 RZI = (STI+STI)*RAZ 179 BRY(2) = 1.0D0/BRY(1) 180 BRY(3) = D1MACH(2) 181 S1R = CYR(1) 182 S1I = CYI(1) 183 S2R = CYR(2) 184 S2I = CYI(2) 185 C1R = CSRR(IFLAG) 186 ASCLE = BRY(IFLAG) 187 K = ND - 2 188 FN = DBLE(FLOAT(K)) 189 DO 100 I=3,ND 190 C2R = S2R 191 C2I = S2I 192 S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I) 193 S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R) 194 S1R = C2R 195 S1I = C2I 196 C2R = S2R*C1R 197 C2I = S2I*C1R 198 YR(K) = C2R 199 YI(K) = C2I 200 K = K - 1 201 FN = FN - 1.0D0 202 IF (IFLAG.GE.3) GO TO 100 203 STR = DABS(C2R) 204 STI = DABS(C2I) 205 C2M = DMAX1(STR,STI) 206 IF (C2M.LE.ASCLE) GO TO 100 207 IFLAG = IFLAG + 1 208 ASCLE = BRY(IFLAG) 209 S1R = S1R*C1R 210 S1I = S1I*C1R 211 S2R = C2R 212 S2I = C2I 213 S1R = S1R*CSSR(IFLAG) 214 S1I = S1I*CSSR(IFLAG) 215 S2R = S2R*CSSR(IFLAG) 216 S2I = S2I*CSSR(IFLAG) 217 C1R = CSRR(IFLAG) 218 100 CONTINUE 219 110 CONTINUE 220 RETURN 221 120 CONTINUE 222 IF (RS1.GT.0.0D0) GO TO 140 223 C----------------------------------------------------------------------- 224 C SET UNDERFLOW AND UPDATE PARAMETERS 225 C----------------------------------------------------------------------- 226 YR(ND) = ZEROR 227 YI(ND) = ZEROI 228 NZ = NZ + 1 229 ND = ND - 1 230 IF (ND.EQ.0) GO TO 110 231 CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM) 232 IF (NUF.LT.0) GO TO 140 233 ND = ND - NUF 234 NZ = NZ + NUF 235 IF (ND.EQ.0) GO TO 110 236 FN = FNU + DBLE(FLOAT(ND-1)) 237 IF (FN.LT.FNUL) GO TO 130 238 C FN = CIDI 239 C J = NUF + 1 240 C K = MOD(J,4) + 1 241 C S1R = CIPR(K) 242 C S1I = CIPI(K) 243 C IF (FN.LT.0.0D0) S1I = -S1I 244 C STR = C2R*S1R - C2I*S1I 245 C C2I = C2R*S1I + C2I*S1R 246 C C2R = STR 247 IN = INU + ND - 1 248 IN = MOD(IN,4) + 1 249 C2R = CAR*CIPR(IN) - SAR*CIPI(IN) 250 C2I = CAR*CIPI(IN) + SAR*CIPR(IN) 251 IF (ZI.LE.0.0D0) C2I = -C2I 252 GO TO 40 253 130 CONTINUE 254 NLAST = ND 255 RETURN 256 140 CONTINUE 257 NZ = -1 258 RETURN 259 150 CONTINUE 260 IF (RS1.GT.0.0D0) GO TO 140 261 NZ = N 262 DO 160 I=1,N 263 YR(I) = ZEROR 264 YI(I) = ZEROI 265 160 CONTINUE 266 RETURN 267 END