github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/internal/amos/amoslib/zbesi.f (about) 1 SUBROUTINE ZBESI(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR) 2 C***BEGIN PROLOGUE ZBESI 3 C***DATE WRITTEN 830501 (YYMMDD) 4 C***REVISION DATE 890801 (YYMMDD) 5 C***CATEGORY NO. B5K 6 C***KEYWORDS I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION, 7 C MODIFIED BESSEL FUNCTION OF THE FIRST KIND 8 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES 9 C***PURPOSE TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT 10 C***DESCRIPTION 11 C 12 C ***A DOUBLE PRECISION ROUTINE*** 13 C ON KODE=1, ZBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX 14 C BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE 15 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE 16 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, ZBESI RETURNS THE SCALED 17 C FUNCTIONS 18 C 19 C CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z) J = 1,...,N , X=REAL(Z) 20 C 21 C WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND 22 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION 23 C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS 24 C (REF. 1). 25 C 26 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION 27 C ZR,ZI - Z=CMPLX(ZR,ZI), -PI.LT.ARG(Z).LE.PI 28 C FNU - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0D0 29 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION 30 C KODE= 1 RETURNS 31 C CY(J)=I(FNU+J-1,Z), J=1,...,N 32 C = 2 RETURNS 33 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N 34 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1 35 C 36 C OUTPUT CYR,CYI ARE DOUBLE PRECISION 37 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS 38 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE 39 C CY(J)=I(FNU+J-1,Z) OR 40 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)) J=1,...,N 41 C DEPENDING ON KODE, X=REAL(Z) 42 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW, 43 C NZ= 0 , NORMAL RETURN 44 C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO 45 C TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0) 46 C J = N-NZ+1,...,N 47 C IERR - ERROR FLAG 48 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED 49 C IERR=1, INPUT ERROR - NO COMPUTATION 50 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) TOO 51 C LARGE ON KODE=1 52 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE 53 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT 54 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE 55 C ACCURACY 56 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA- 57 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI- 58 C CANCE BY ARGUMENT REDUCTION 59 C IERR=5, ERROR - NO COMPUTATION, 60 C ALGORITHM TERMINATION CONDITION NOT MET 61 C 62 C***LONG DESCRIPTION 63 C 64 C THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR 65 C SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z), 66 C THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A 67 C NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE 68 C UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z) 69 C FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE 70 C SEQUENCES OR REDUCE ORDERS WHEN NECESSARY. 71 C 72 C THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND 73 C CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA 74 C 75 C I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z) REAL(Z).GT.0.0 76 C M = +I OR -I, I**2=-1 77 C 78 C FOR NEGATIVE ORDERS,THE FORMULA 79 C 80 C I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z) 81 C 82 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE 83 C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE 84 C INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE 85 C NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER, 86 C K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF 87 C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY 88 C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN 89 C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE, 90 C LARGE MEANS FNU.GT.CABS(Z). 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 ZBINU,I1MACH,D1MACH 151 C***END PROLOGUE ZBESI 152 C COMPLEX CONE,CSGN,CW,CY,CZERO,Z,ZN 153 DOUBLE PRECISION AA, ALIM, ARG, CONEI, CONER, CSGNI, CSGNR, CYI, 154 * CYR, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR, 155 * ZR, D1MACH, AZ, BB, FN, ZABS, ASCLE, RTOL, ATOL, STI 156 INTEGER I, IERR, INU, K, KODE, K1,K2,N,NZ,NN, I1MACH 157 DIMENSION CYR(N), CYI(N) 158 DATA PI /3.14159265358979324D0/ 159 DATA CONER, CONEI /1.0D0,0.0D0/ 160 C 161 C***FIRST EXECUTABLE STATEMENT ZBESI 162 IERR = 0 163 NZ=0 164 IF (FNU.LT.0.0D0) IERR=1 165 IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 166 IF (N.LT.1) IERR=1 167 IF (IERR.NE.0) RETURN 168 C----------------------------------------------------------------------- 169 C SET PARAMETERS RELATED TO MACHINE CONSTANTS. 170 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. 171 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. 172 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND 173 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR 174 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. 175 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z. 176 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG). 177 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU. 178 C----------------------------------------------------------------------- 179 TOL = DMAX1(D1MACH(4),1.0D-18) 180 K1 = I1MACH(15) 181 K2 = I1MACH(16) 182 R1M5 = D1MACH(5) 183 K = MIN0(IABS(K1),IABS(K2)) 184 ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0) 185 K1 = I1MACH(14) - 1 186 AA = R1M5*DBLE(FLOAT(K1)) 187 DIG = DMIN1(AA,18.0D0) 188 AA = AA*2.303D0 189 ALIM = ELIM + DMAX1(-AA,-41.45D0) 190 RL = 1.2D0*DIG + 3.0D0 191 FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0) 192 C----------------------------------------------------------------------------- 193 C TEST FOR PROPER RANGE 194 C----------------------------------------------------------------------- 195 AZ = ZABS(CMPLX(ZR,ZI,kind=KIND(1.0D0))) 196 FN = FNU+DBLE(FLOAT(N-1)) 197 AA = 0.5D0/TOL 198 BB=DBLE(FLOAT(I1MACH(9)))*0.5D0 199 AA = DMIN1(AA,BB) 200 IF (AZ.GT.AA) GO TO 260 201 IF (FN.GT.AA) GO TO 260 202 AA = DSQRT(AA) 203 IF (AZ.GT.AA) IERR=3 204 IF (FN.GT.AA) IERR=3 205 ZNR = ZR 206 ZNI = ZI 207 CSGNR = CONER 208 CSGNI = CONEI 209 IF (ZR.GE.0.0D0) GO TO 40 210 ZNR = -ZR 211 ZNI = -ZI 212 C----------------------------------------------------------------------- 213 C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE 214 C WHEN FNU IS LARGE 215 C----------------------------------------------------------------------- 216 INU = INT(SNGL(FNU)) 217 ARG = (FNU-DBLE(FLOAT(INU)))*PI 218 IF (ZI.LT.0.0D0) ARG = -ARG 219 CSGNR = DCOS(ARG) 220 CSGNI = DSIN(ARG) 221 IF (MOD(INU,2).EQ.0) GO TO 40 222 CSGNR = -CSGNR 223 CSGNI = -CSGNI 224 40 CONTINUE 225 CALL ZBINU(ZNR, ZNI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, 226 * ELIM, ALIM) 227 IF (NZ.LT.0) GO TO 120 228 IF (ZR.GE.0.0D0) RETURN 229 C----------------------------------------------------------------------- 230 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE 231 C----------------------------------------------------------------------- 232 NN = N - NZ 233 IF (NN.EQ.0) RETURN 234 RTOL = 1.0D0/TOL 235 ASCLE = D1MACH(1)*RTOL*1.0D+3 236 DO 50 I=1,NN 237 C STR = CYR(I)*CSGNR - CYI(I)*CSGNI 238 C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR 239 C CYR(I) = STR 240 AA = CYR(I) 241 BB = CYI(I) 242 ATOL = 1.0D0 243 IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 55 244 AA = AA*RTOL 245 BB = BB*RTOL 246 ATOL = TOL 247 55 CONTINUE 248 STR = AA*CSGNR - BB*CSGNI 249 STI = AA*CSGNI + BB*CSGNR 250 CYR(I) = STR*ATOL 251 CYI(I) = STI*ATOL 252 CSGNR = -CSGNR 253 CSGNI = -CSGNI 254 50 CONTINUE 255 RETURN 256 120 CONTINUE 257 IF(NZ.EQ.(-2)) GO TO 130 258 NZ = 0 259 IERR=2 260 RETURN 261 130 CONTINUE 262 NZ=0 263 IERR=5 264 RETURN 265 260 CONTINUE 266 NZ=0 267 IERR=4 268 RETURN 269 END