github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zbesk.f (about) 1 SUBROUTINE ZBESK(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) 2 C***BEGIN PROLOGUE ZBESK 3 C***DATE WRITTEN 830501 (YYMMDD) 4 C***REVISION DATE 890801 (YYMMDD) 5 C***CATEGORY NO. B5K 6 C***KEYWORDS K-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION, 7 C MODIFIED BESSEL FUNCTION OF THE SECOND KIND, 8 C BESSEL FUNCTION OF THE THIRD KIND 9 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 10 C***PURPOSE TO COMPUTE K-BESSEL FUNCTIONS OF COMPLEX ARGUMENT 11 C***DESCRIPTION 12 C 13 C ***A DOUBLE PRECISION ROUTINE*** 14 C 15 C ON KODE=1, CBESK COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 16 C BESSEL FUNCTIONS CY(J)=K(FNU+J-1,Z) FOR REAL, NONNEGATIVE 17 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z.NE.CMPLX(0.0,0.0) 18 C IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESK 19 C RETURNS THE SCALED K FUNCTIONS, 20 C 21 C CY(J)=EXP(Z)*K(FNU+J-1,Z) , J=1,...,N, 22 C 23 C WHICH REMOVE THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND 24 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND 25 C NOTATION ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL 26 C FUNCTIONS (REF. 1). 27 C 28 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 29 C ZR,ZI - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0), 30 C -PI.LT.ARG(Z).LE.PI 31 C FNU - ORDER OF INITIAL K FUNCTION, FNU.GE.0.0D0 32 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 33 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 34 C KODE= 1 RETURNS 35 C CY(I)=K(FNU+I-1,Z), I=1,...,N 36 C = 2 RETURNS 37 C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N 38 C 39 C OUTPUT CYR,CYI ARE DOUBLE PRECISION 40 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 41 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 42 C CY(I)=K(FNU+I-1,Z), I=1,...,N OR 43 C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N 44 C DEPENDING ON KODE 45 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW. 46 C NZ= 0 , NORMAL RETURN 47 C NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO DUE 48 C TO UNDERFLOW, CY(I)=CMPLX(0.0D0,0.0D0), 49 C I=1,...,N WHEN X.GE.0.0. WHEN X.LT.0.0 50 C NZ STATES ONLY THE NUMBER OF UNDERFLOWS 51 C IN THE SEQUENCE. 52 C 53 C IERR - ERROR FLAG 54 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 55 C IERR=1, INPUT ERROR - NO COMPUTATION 56 C IERR=2, OVERFLOW - NO COMPUTATION, FNU IS 57 C TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH 58 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 59 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 60 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 61 C ACCURACY 62 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 63 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 64 C CANCE BY ARGUMENT REDUCTION 65 C IERR=5, ERROR - NO COMPUTATION, 66 C ALGORITHM TERMINATION CONDITION NOT MET 67 C 68 C***LONG DESCRIPTION 69 C 70 C EQUATIONS OF THE REFERENCE ARE IMPLEMENTED FOR SMALL ORDERS 71 C DNU AND DNU+1.0 IN THE RIGHT HALF PLANE X.GE.0.0. FORWARD 72 C RECURRENCE GENERATES HIGHER ORDERS. K IS CONTINUED TO THE LEFT 73 C HALF PLANE BY THE RELATION 74 C 75 C K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z) 76 C MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1 77 C 78 C WHERE I(FNU,Z) IS THE I BESSEL FUNCTION. 79 C 80 C FOR LARGE ORDERS, FNU.GT.FNUL, THE K FUNCTION IS COMPUTED 81 C BY MEANS OF ITS UNIFORM ASYMPTOTIC EXPANSIONS. 82 C 83 C FOR NEGATIVE ORDERS, THE FORMULA 84 C 85 C K(-FNU,Z) = K(FNU,Z) 86 C 87 C CAN BE USED. 88 C 89 C CBESK ASSUMES THAT A SIGNIFICANT DIGIT SINH(X) FUNCTION IS 90 C AVAILABLE. 91 C 92 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE- 93 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS 94 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. 95 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN 96 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG 97 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS 98 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION. 99 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS 100 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS 101 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE 102 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS 103 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3 104 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION 105 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION 106 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN 107 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT 108 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS 109 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC. 110 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES. 111 C 112 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX 113 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT 114 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE- 115 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE 116 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))), 117 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF 118 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY 119 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN 120 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY 121 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER 122 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K, 123 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS 124 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER 125 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY 126 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER 127 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE 128 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES, 129 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P, 130 C OR -PI/2+P. 131 C 132 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ 133 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF 134 C COMMERCE, 1955. 135 C 136 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 137 C BY D. E. AMOS, SAND83-0083, MAY, 1983. 138 C 139 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT 140 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983. 141 C 142 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 143 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85- 144 C 1018, MAY, 1985 145 C 146 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX 147 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS. 148 C MATH. SOFTWARE, 1986 149 C 150 C***ROUTINES CALLED ZACON,ZBKNU,ZBUNK,ZUOIK,ZABS,I1MACH,D1MACH 151 C***END PROLOGUE ZBESK 152 C 153 C COMPLEX CY,Z 154 DOUBLE PRECISION AA, ALIM, ALN, ARG, AZ, CYI, CYR, DIG, ELIM, FN, 155 * FNU, FNUL, RL, R1M5, TOL, UFL, ZI, ZR, D1MACH, ZABS, BB 156 INTEGER IERR, K, KODE, K1, K2, MR, N, NN, NUF, NW, NZ, I1MACH 157 DIMENSION CYR(N), CYI(N) 158 C***FIRST EXECUTABLE STATEMENT ZBESK 159 IERR = 0 160 NZ=0 161 IF (ZI.EQ.0.0E0 .AND. ZR.EQ.0.0E0) IERR=1 162 IF (FNU.LT.0.0D0) IERR=1 163 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 164 IF (N.LT.1) IERR=1 165 IF (IERR.NE.0) RETURN 166 NN = N 167 C----------------------------------------------------------------------- 168 C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 169 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 170 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 171 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 172 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 173 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 174 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 175 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 176 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU 177 C----------------------------------------------------------------------- 178 TOL = DMAX1(D1MACH(4),1.0D-18) 179 K1 = I1MACH(15) 180 K2 = I1MACH(16) 181 R1M5 = D1MACH(5) 182 K = MIN0(IABS(K1),IABS(K2)) 183 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 184 K1 = I1MACH(14) - 1 185 AA = R1M5*DBLE(FLOAT(K1)) 186 DIG = DMIN1(AA,18.0D0) 187 AA = AA*2.303D0 188 ALIM = ELIM + DMAX1(-AA,-41.45D0) 189 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 190 RL = 1.2D0*DIG + 3.0D0 191 C----------------------------------------------------------------------------- 192 C TEST FOR PROPER RANGE 193 C----------------------------------------------------------------------- 194 AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 195 FN = FNU + DBLE(FLOAT(NN-1)) 196 AA = 0.5D0/TOL 197 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 198 AA = DMIN1(AA,BB) 199 IF (AZ.GT.AA) GO TO 260 200 IF (FN.GT.AA) GO TO 260 201 AA = DSQRT(AA) 202 IF (AZ.GT.AA) IERR=3 203 IF (FN.GT.AA) IERR=3 204 C----------------------------------------------------------------------- 205 C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE 206 C----------------------------------------------------------------------- 207 C UFL = DEXP(-ELIM) 208 UFL = D1MACH(1)*1.0D+3 209 IF (AZ.LT.UFL) GO TO 180 210 IF (FNU.GT.FNUL) GO TO 80 211 IF (FN.LE.1.0D0) GO TO 60 212 IF (FN.GT.2.0D0) GO TO 50 213 IF (AZ.GT.TOL) GO TO 60 214 ARG = 0.5D0*AZ 215 ALN = -FN*DLOG(ARG) 216 IF (ALN.GT.ELIM) GO TO 180 217 GO TO 60 218 50 CONTINUE 219 CALL ZUOIK(ZR, ZI, FNU, KODE, 2, NN, CYR, CYI, NUF, TOL, ELIM, 220 * ALIM) 221 IF (NUF.LT.0) GO TO 180 222 NZ = NZ + NUF 223 NN = NN - NUF 224 C----------------------------------------------------------------------- 225 C HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK 226 C IF NUF=NN, THEN CY(I)=CZERO FOR ALL I 227 C----------------------------------------------------------------------- 228 IF (NN.EQ.0) GO TO 100 229 60 CONTINUE 230 IF (ZR.LT.0.0D0) GO TO 70 231 C----------------------------------------------------------------------- 232 C RIGHT HALF PLANE COMPUTATION, REAL(Z).GE.0. 233 C----------------------------------------------------------------------- 234 CALL ZBKNU(ZR, ZI, FNU, KODE, NN, CYR, CYI, NW, TOL, ELIM, ALIM) 235 IF (NW.LT.0) GO TO 200 236 NZ=NW 237 RETURN 238 C----------------------------------------------------------------------- 239 C LEFT HALF PLANE COMPUTATION 240 C PI/2.LT.ARG(Z).LE.PI AND -PI.LT.ARG(Z).LT.-PI/2. 241 C----------------------------------------------------------------------- 242 70 CONTINUE 243 IF (NZ.NE.0) GO TO 180 244 MR = 1 245 IF (ZI.LT.0.0D0) MR = -1 246 CALL ZACON(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, NW, RL, FNUL, 247 * TOL, ELIM, ALIM) 248 IF (NW.LT.0) GO TO 200 249 NZ=NW 250 RETURN 251 C----------------------------------------------------------------------- 252 C UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL 253 C----------------------------------------------------------------------- 254 80 CONTINUE 255 MR = 0 256 IF (ZR.GE.0.0D0) GO TO 90 257 MR = 1 258 IF (ZI.LT.0.0D0) MR = -1 259 90 CONTINUE 260 CALL ZBUNK(ZR, ZI, FNU, KODE, MR, NN, CYR, CYI, NW, TOL, ELIM, 261 * ALIM) 262 IF (NW.LT.0) GO TO 200 263 NZ = NZ + NW 264 RETURN 265 100 CONTINUE 266 IF (ZR.LT.0.0D0) GO TO 180 267 RETURN 268 180 CONTINUE 269 NZ = 0 270 IERR=2 271 RETURN 272 200 CONTINUE 273 IF(NW.EQ.(-1)) GO TO 180 274 NZ=0 275 IERR=5 276 RETURN 277 260 CONTINUE 278 NZ=0 279 IERR=4 280 RETURN 281 END