! Licensed under the ACM Software License Agreement ! Copyright © 1970–2012 Association for Computing Machinery (ACM) ! See ColSpecF LICENSE file for details. MODULE calgo_683 !* # CALGO 683 ! Algorithm 683. ! ! Procedures: ! ! - `cexint`: Exponential integral \(\mathrm{E}_n(z)\) ! ! Private procedures: ! ! - `psixn`: Digamma function \(\psi(n)\) for a positive integer \(n\) ! ! ## Author ! Donald E. Amos ! ! ## History ! - 1987-05-15 - Donald E. Amos ! - Original code. ! - 1999-12-28 - Donald E. Amos ! - F77 code distributed by ACM: <https://calgo.acm.org/683.zip> ! - YYYY-mm-dd - Alan Miller ! - F90 code adaptation by Alan Miller: ! <https://jblevins.org/mirror/amiller/toms683.f90> ! - 2025-06-05 - Rodrigo Castro (GitHub: rodpcastro) ! - Removed procedures `cexqad` and `fqcex`, which were originally used for ! testing. ! - 2025-06-06 - Rodrigo Castro (GitHub: rodpcastro) ! - Created abstract interface for single-variable function, which is used by `g8` ! and `gaus8`. ! - Replaced `dp` (double precision) by `wp` (working precision). ! - Fixed typo at `gaus8`: ! - `anib = LOG10(DBLE(RADIX(0.0_wp))) * k / 0.30102000_wp` ! - to `anib = LOG10(DBLE(RADIX(0.0_wp))) * k / 0.30103000_wp`. ! - 2025-07-27 - Rodrigo Castro (GitHub: rodpcastro) ! - Removed procedures `g8` and `gaus8`, which were originally used for testing. ! - Turned `psixn` into a private function. ! ! ## References ! 1. Donald E. Amos. 1990. Algorithms 683: a portable FORTRAN subroutine for ! exponential integrals of a complex argument. ACM Trans. Math. Softw. 16, !* 2 (June 1990), 178–182. <https://doi.org/10.1145/78928.78934> use csf_kinds, only: wp IMPLICIT NONE PRIVATE PUBLIC :: cexint REAL(wp), SAVE :: fnm, gln, x, y INTEGER, SAVE :: iflag CONTAINS SUBROUTINE cexint(z, n, kode, tol, m, cy, ierr) !! CALGO 683 Exponential integral \(\mathrm{E}_n(z)\). ! !! \(n \geq 1,\thinspace !! \lbrace z \in \mathbb{C} \mid -\pi \lt \arg(z) \leq \pi \rbrace \) ! !! The description of all parameters and outputs can be found in the source code !! docstring. ! ! ON KODE=1, CEXINT COMPUTES AN M MEMBER SEQUENCE OF COMPLEX(wp) ! EXPONENTIAL INTEGRALS CY(J)=E(N+J-1,Z), J=1,...,M, FOR ! POSITIVE ORDERS N,...,N+M-1 AND COMPLEX(wp) Z IN THE CUT PLANE ! -PI < ARG(Z) <= PI (N=1 AND Z=CMPLX(0.0,0.0) CANNOT HOLD AT ! THE SAME TIME). ON KODE=2, CEXINT COMPUTES SCALED FUNCTIONS ! ! CY(J)=E(N+J-1,Z)*CEXP(Z), J=1,...,M, ! ! WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE LEFT AND ! RIGHT HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN THE ! NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1). ! ! INPUT ! Z - Z=CMPLX(X,Y), -PI < ARG(Z) <= PI ! N - INTEGER ORDER OF INITIAL E FUNCTION, N=1,2,... ! (N=1 AND Z=CMPLX(0.0,0.0) IS AN ERROR) ! KODE - A PARAMETER TO INDICATE THE SCALING OPTION ! KODE= 1 RETURNS ! CY(J)=E(N+J-1,Z), J=1,...,M ! = 2 RETURNS ! CY(J)=E(N+J-1,Z)*CEXP(Z), J=1,...,M ! TOL - PRECISION (ACCURACY) DESIRED FOR THE SEQUENCE, ! URND <= TOL <= 1.0E-3, WHERE URND IS LIMITED BY ! URND = MAX(UNIT ROUNDOFF,1.0E-18) AND UNIT ! ROUNDOFF = EPSILON(0.0_wp) ! M - NUMBER OF E FUNCTIONS IN THE SEQUENCE, M >= 1 ! ! OUTPUT ! CY - A COMPLEX(wp) VECTOR WHOSE FIRST M COMPONENTS CONTAIN ! VALUES FOR THE SEQUENCE ! CY(J)=E(N+J-1,Z) OR ! CY(J)=E(N+J-1,Z)*CEXP(Z), J=1,...,M ! DEPENDING ON KODE. ! IERR - ERROR FLAG ! IERR=0, NORMAL RETURN - COMPUTATION COMPLETED ! IERR=1, INPUT ERROR - NO COMPUTATION ! IERR=2, UNDERFLOW - FIRST M COMPONENTS OF CY ! SET TO ZERO, CY(J)=CMPLX(0.0,0.0), J=1,M, ! REAL(Z) > 0.0 TOO LARGE ON KODE=1 ! IERR=3, OVERFLOW - NO COMPUTATION, ! REAL(Z) < 0.0 TOO SMALL ON KODE=1 ! IERR=4, CABS(Z) OR N+M-1 LARGE - COMPUTATION DONE ! BUT LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION ! PRODUCE LESS THAN HALF OF MACHINE ACCURACY ! IERR=5, CABS(Z) OR N+M-1 TOO LARGE - NO COMPUTATION BECAUSE ! OF COMPLETE LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION ! IERR=6, ERROR - NO COMPUTATION, ! ALGORITHM TERMINATION CONDITION NOT MET. ! SEE LONG DESCRIPTION ABOUT PARAMETER ICMAX. ! IERR=7, ERROR - NO COMPUTATION, ! DISCRIMINATION ERROR. THIS CONDITION SHOULD NEVER OCCUR. ! ! LONG DESCRIPTION ! ! CEXINT USES A COMBINATION OF POWER SERIES AND BACKWARD RECURRENCE ! DESCRIBED IN REF. 2 FOR THE COMPLEX(wp) Z PLANE EXCEPT FOR A STRIP ! 2*YB WIDE ABOUT THE NEGATIVE REAL AXIS, WHERE ANALYTIC CONTINUATION IS ! CARRIED OUT BY LIMITED USE OF TAYLOR SERIES. ! THE SWITCH FROM BACKWARD RECURRENCE TO TAYLOR SERIES IS NECESSARY BECAUSE ! BACKWARD RECURRENCE IS SLOWLY CONVERGENT NEAR THE NEGATIVE REAL AXIS. ! THE BOUNDARIES Y=-YB AND Y=YB WERE DETERMINED SO THAT BACKWARD RECURRENCE ! WOULD CONVERGE EASILY WITH N AS LARGE AS 100 AND TOL AS SMALL AS 1.0e-18. ! SUBROUTINE CEXENZ DOES THE BACKWARD RECURRENCE AND SUBROUTINE CACEXI DOES ! THE ANALYTIC CONTINUATION. TO START THE CONTINUATION, CACEXI CALLS CEXENZ ! WITH ZB=CMPLX(X,YB). ! IF CEXENZ RETURNS IERR=6, THEN YB IS INCREASED BY 0.5 UNTIL CEXENZ RETURNS ! IERR=0 OR 10 TRIES, WHICHEVER COMES FIRST. ! WHEN IERR=0, THEN THE ANALYTIC CONTINUATION PROCEEDS VERTICALLY DOWN FROM ! ZB=CMPLX(X,YB) TO Z=CMPLX(X,Y), 0 <= Y < YB. ! CONJUGATION IS USED FOR Y < 0. YB INCREASES AS TOL DECREASES TO KEEP ! CONVERGENCE RATES UP AND RECURRENCE DOWN. ! ! PARAMETER ICDIM=250 ALLOCATES STORAGE FOR THE COEFFICIENTS OF THE BACKWARD ! RECURRENCE ALGORITHM. IF THE ALGORITHM TERMINATION CONDITION IS NOT MET ! IN ICDIM STEPS, THEN RECURRENCE PROCEEDS WITH NO ADDITIONAL STORAGE UNTIL ! THE TERMINATION CONDITION IS MET OR THE LIMIT ICMAX=2000 IS EXCEEDED. ! THE PURPOSE OF STORAGE IS TO MAKE THE ALGORITHM MORE EFFICIENT. ! THE TERMINATION CONDITION IS MET IN LESS THAN 250 STEPS OVER MOST OF ! THE COMPLEX(wp) PLANE EXCLUDING THE STRIP ABS(Y) < YB, X < 0. ! EXCEPTIONS TO THIS RULE ARE GENERATED NEAR STRIP BOUNDARIES WHEN N+M-1 ! AND ABS(Z) ARE LARGE AND NEARLY EQUAL. IN THESE CASES, THE CONVERGENCE ! IS VERY SLOW AND ADDITIONAL RECURRENCE (UP TO ICMAX) MUST BE USED. ! ON THE OTHERHAND, THESE REGIONS OF SLOW CONVERGENCE ARE KEPT SMALL BY ! ADJUSTING YB AS A FUNCTION OF TOL. THESE REGIONS COULD BE ELIMINATED ! ENTIRELY BY MAKING YB SUFFICIENTLY LARGE, BUT THE EXPENSE AND INSTABILITY ! OF CONTINUATION BY TAYLOR SERIES NOT ONLY GOES UP, BUT THE COMPUTATIONAL ! EXPENSE BECOMES EXCESSIVELY LARGE IN OTHER PARTS OF THE LEFT HALF PLANE ! (Y < YB) WHERE THE BACKWARD RECURRENCE ALGORITHM WOULD CONVERGE RAPIDLY. ! ! DERIVATIVES FOR SUCCESSIVE POWER SERIES ARE NOT COMPUTED BY EVALUATING ! DERIVATIVES OF A PREVIOUS POWER SERIES. BECAUSE OF THE RELATION ! ! (1) DE(N,Z)/DZ = - E(N-1,Z), ! ! SUCCESSIVE DERIVATIVES AT Z ARE GIVEN BY LOWER ORDER FUNCTIONS AND CAN BE ! COMPUTED IN A STABLE FASHION BY BACKWARD RECURRENCE USING (2) PROVIDED ! THAT THE BEGINNING ORDER NUB IS SMALLER THAN THE ARGUMENT. ! TO ACHIEVE THIS FOR ALL INTERMEDIATE VALUES ZZ BETWEEN ZB AND Z, WE TAKE ! NUB=MINO(N+M-1,INT(CABS(Z)+0.5)). ! TO START, E(NUB,ZB) IS EVALUATED BY THE BACKWARD RECURRENCE ALGORITHM OF ! REF. 3. TO CONTINUE THE FUNCTION FROM ZB TO Z VIA INTERMEDIATE VALUES ZZ, ! DERIVATIVES OF E(NUB,ZB) ARE COMPUTED BY BACKWARD RECURRENCE ON (2). ! THIS ALLOWS A STEP (NO LARGER THAN 0.5) TO ZZ FOR E(NUB,ZZ) USING THE ! TAYLOR SERIES ABOUT ZB. NOW, WE APPLY (2) AGAIN STARTING AT E(NUB,ZZ) FOR ! THE DERIVATIVES AT ZZ AND TAKE ANOTHER STEP, ETC. NOTICE THAT THE ! STABILITY CONDITION FOR BACKWARD RECURRENCE, NUB <= ABS(Z) <= ABS(ZZ) ! <= ABS(ZB), IS SATISFIED FOR ALL INTERMEDIATE VALUES ZZ. THE FINAL ! SEQUENCE FOR ORDERS N,...,N+M-1 IS GENERATED FROM (2) BY BACKWARD ! RECURRENCE, FORWARD RECURRENCE OR BOTH ONCE E(NUB,Z) HAS BEEN COMPUTED. ! ! RECURRENCE WITH THE RELATION ! ! (2) N*E(N+1,Z) + Z*E(N,Z) = CEXP(-Z) ! ! IN A DIRECTION AWAY FROM THE INTEGER CLOSEST TO ABS(Z) IS STABLE. ! FOR NEGATIVE ORDERS, THE RECURRENCE ! ! E( 0,Z) = CEXP(-Z)/Z ! E(-N,Z) = ( CEXP(-Z)+N*E(-N+1,Z) )/Z ,N=1,2,... ! ! IS NUMERICALLY STABLE FOR ALL Z. ! ! THE (CAPITAL) SINE AND COSINE INTEGRALS CAN BE COMPUTED FROM ! ! SI(Z) = (E(1,I*Z)-E(1,-I*Z))/(2*I) + PI/2 ! CI(Z) = -(E(1,I*Z)+E(1,-I*Z))/2 ! ! IN -PI/2 < ARG(Z) <= PI/2, (I**2=-1), WHILE THE PRINCIPAL ! VALUED EXPONENTIAL INTEGRAL EI(X) CAN BE COMPUTED FROM ! ! EI( X) = -(E(1,-X+I*0)+E(1,-X-I*0))/2 = -REAL(E(1,-X)) ! EI(-X) = -REAL(E(1,X)) ! ! FOR X > 0.0 TO AN ACCURACY TOL. IF Z = X > 0 THEN THE REAL SINE AND ! COSINE INTEGRALS ARE GIVEN BY ! ! SI(X) = AIMAG(E(1,I*X)) + PI/2 ! CI(X) = -REAL(E(1,I*X)) . ! ! THE ANALYTIC CONTINUATION TO OTHER SHEETS CAN BE DONE BY THE RELATIONS ! ! E(N,Z*CEXP(2*PI*M*I)) = E(N,Z) - 2*PI*M*I*(-Z)**(N-1)/(N-1)! ! ! WHERE M=+1 OR M=-1 AND I**2=-1. ! ! REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ AND ! I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF COMMERCE, 1955. ! ! COMPUTATION OF EXPONENTIAL INTEGRALS OF COMPLEX(wp) ARGUMENT ! BY D. E. AMOS, ACM TRANS. MATH. SOFTWARE ! ! COMPUTATION OF EXPONENTIAL INTEGRALS BY D. E. AMOS, ACM TRANS. ! MATH. SOFTWARE, VOL 6, NO. 3 SEPTEMBER 1980, PP. 365-377; ! ALGORITHM 556, EXPONENTIAL INTEGRALS, PP. 420-428. ! ! REMARK ON ALGORITHM 556 ! BY D. E. AMOS, ACM TRANS. MATH. SOFTWARE, VOL 9, NO. 4 ! DECEMBER 1983, P. 525. ! ! UNIFORM ASYMPTOTIC EXPANSIONS FOR EXPONENTIAL INTEGRALS E(N,X) ! AND BICKLEY FUNCTIONS KI(N,X) BY D. E. AMOS, ACM TRANS. MATH. ! SOFTWARE, VOL 9, NO. 4, DECEMBER. 1983, PP. 467-479; ! ALGORITHM 609, A PORTABLE FORTRAN SUBROUTINE FOR BICKLEY ! FUNCTIONS KI(N,X), PP. 480-493. ! ! ROUTINES CALLED CACEXI, CEXENZ COMPLEX(wp), INTENT(IN) :: z INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: kode REAL(wp), INTENT(IN) :: tol INTEGER, INTENT(IN) :: m COMPLEX(wp), INTENT(OUT) :: cy(m) INTEGER, INTENT(OUT) :: ierr INTEGER :: i, k, k1, k2 REAL(wp) :: aa, alim, az, bb, ca(250), d, elim, & fn, rbry, r1m5, urnd, x, y, yb, htol COMPLEX(wp) :: cb(250) !----------------------------------------------------------------------- ! DIMENSION CA(ICDIM),CB(ICDIM) !----------------------------------------------------------------------- INTEGER, PARAMETER :: icdim = 250 ! FIRST EXECUTABLE STATEMENT CEXINT ierr = 0 x = REAL(z, KIND=wp) y = AIMAG(z) IF (x == 0.0_wp .AND. y == 0.0_wp .AND. n == 1) ierr = 1 IF (n < 1) ierr = 1 IF (kode < 1 .OR. kode > 2) ierr = 1 IF (m < 1) ierr = 1 urnd = MAX(EPSILON(0.0_wp), 1.0E-18_wp) IF (tol < urnd .OR. tol > 1.0E-3_wp) ierr = 1 IF (ierr /= 0) RETURN IF (x /= 0.0_wp .OR. y /= 0.0_wp .OR. n <= 1) THEN !----------------------------------------------------------------------- ! SET PARAMETERS RELATED TO MACHINE CONSTANTS. ! URND IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18. ! ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT. ! EXP(-ELIM) < EXP(-ALIM)=EXP(-ELIM)/URND AND ! EXP(ELIM) > EXP(ALIM)=EXP(ELIM)*URND ARE INTERVALS NEAR ! UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE. !----------------------------------------------------------------------- k1 = MINEXPONENT(0.0_wp) k2 = MAXEXPONENT(0.0_wp) r1m5 = LOG10(DBLE(RADIX(0.0_wp))) k = MIN(ABS(k1),ABS(k2)) elim = 2.303_wp * (k*r1m5 - 3.0_wp) k1 = DIGITS(0.0_wp) - 1 aa = r1m5 * k1 aa = aa * 2.303 alim = elim + MAX(-aa, -41.45_wp) rbry = 2.0 IF (urnd > 1.0E-8_wp) rbry = 1.0_wp !----------------------------------------------------------------------- ! TEST VARIABLES FOR RANGE. ABS(Z) CANNOT BE LARGER THAN THE ARGUMENT ! OF THE INT( ) FUNCTION. !----------------------------------------------------------------------- az = ABS(z) fn = n+m-1 aa = 0.5_wp / urnd bb = HUGE(0.0_wp) * 0.5_wp aa = MIN(aa,bb) IF (az > aa) GO TO 20 IF (fn > aa) GO TO 20 aa = SQRT(aa) IF (az > aa) ierr = 4 IF (fn > aa) ierr = 4 IF (x >= 0.0_wp) THEN !----------------------------------------------------------------------- ! BACKWARD RECURRENCE FOR THE RIGHT HALF PLANE, X >= 0.0E0 !----------------------------------------------------------------------- CALL cexenz(z, n, kode, m, cy, ierr, rbry, tol, elim, alim, icdim, ca, cb) RETURN END IF IF (az <= rbry) THEN !----------------------------------------------------------------------- ! POWER SERIES FOR ABS(Z) <= RBRY AND X < 0.0E0 !----------------------------------------------------------------------- CALL cexenz(z, n, kode, m, cy, ierr, rbry, tol, elim, alim, icdim, ca, cb) RETURN END IF d = -0.4342945_wp * LOG(tol) yb = 10.5_wp - 0.538460_wp * (18.0_wp-d) IF (ABS(y) >= yb) THEN !----------------------------------------------------------------------- ! BACKWARD RECURRENCE EXTERIOR TO THE STRIP ABS(Y) < YB, X < 0.0 !----------------------------------------------------------------------- htol = 0.125_wp * tol CALL cexenz(z, n, kode, m, cy, ierr, rbry, htol, elim, alim, icdim, ca, cb) RETURN END IF !----------------------------------------------------------------------- ! TAYLOR SERIES IN CACEXI FOR ANALYTIC CONTINUATION !----------------------------------------------------------------------- CALL cacexi(z, n, kode, m, cy, ierr, yb, rbry, tol, elim, alim, icdim, ca, cb ) RETURN END IF DO i = 1, m cy(i) = 1.0_wp / (n+i-2) END DO RETURN 20 ierr = 5 RETURN END SUBROUTINE cexint SUBROUTINE cexenz(z, n, kode, m, cy, ierr, rbry, tol, elim, alim, icdim, ca, cb) ! REFER TO CEXINT ! ! CEXENZ COMPUTES THE EXPONENTIAL INTEGRAL BY MEANS OF POWER SERIES AND A ! BACKWARD RECURRENCE ALGORITHM FOR THE CONFLUENT HYPERGEOMETRIC ! REPRESENTATION. ! ! ROUTINES CALLED PSIXN COMPLEX(wp), INTENT(IN) :: z INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: kode INTEGER, INTENT(IN) :: m COMPLEX(wp), INTENT(OUT) :: cy(m) INTEGER, INTENT(OUT) :: ierr REAL(wp), INTENT(IN) :: rbry REAL(wp), INTENT(IN) :: tol REAL(wp), INTENT(IN) :: elim REAL(wp), INTENT(IN) :: alim INTEGER, INTENT(IN) :: icdim REAL(wp), INTENT(OUT) :: ca(icdim) COMPLEX(wp), INTENT(OUT) :: cb(icdim) INTEGER :: i, ic, icase, ict, ik, ind, iz, jset, & k, kflag, kn, ks, ml, mu, nd, nm REAL(wp) :: aa, aam, aams, aem, ah, ak, ap1, at, az, bk, bt, & dk, ERR, fc, fnm, rtola, tola, x, xtol, y, ck COMPLEX(wp) :: cp1, cp2, cpt, cat, cbt, cy1, cy2, cyy(2), cnorm, & cs, cak, emz, caa, tz, fz, ct, scle, rscle REAL(wp), PARAMETER :: euler = -5.77215664901532861e-1_wp COMPLEX(wp), PARAMETER :: czero = (0.0_wp, 0.0_wp), cone = (1.0_wp, 0.0_wp) INTEGER, PARAMETER :: icmax = 2000 ierr = 0 scle = cone rscle = cone x = REAL(z, KIND=wp) y = AIMAG(z) az = ABS(z) IF (az <= rbry) THEN !----------------------------------------------------------------------- ! SERIES FOR E(N,Z) FOR ABS(Z) <= RBRY !----------------------------------------------------------------------- iz = az + 0.5_wp !----------------------------------------------------------------------- ! ICASE=1 MEANS INTEGER CLOSEST TO ABS(Z) IS 2 AND N=1 ! ICASE=2 MEANS INTEGER CLOSEST TO ABS(Z) IS 0,1, OR 2 AND N >= 2 !----------------------------------------------------------------------- icase = 2 IF (iz > n) icase = 1 nm = n - icase + 1 nd = nm + 1 ind = 3 - icase mu = m - ind ml = 1 ks = nd fnm = nm cs = czero xtol = 0.3333_wp * tol aam = 1.0_wp IF (nd /= 1) THEN aam = 1.0_wp / fnm cs = aam END IF caa = cone aa = 1.0_wp ak = 1.0_wp !----------------------------------------------------------------------- ! LIMIT INDEX I TO IK TO PREVENT UNDERFLOW ON SMALL VALUES OF Z !----------------------------------------------------------------------- ik = 35 IF (az < xtol*aam) ik = 1 DO i = 1, ik at = 1.0_wp / ak caa = -caa * z * at aa = aa * az * at IF (i /= nm) THEN cs = cs - caa / (ak-fnm) ELSE cs = cs + caa * (-LOG(z) + psixn(nd)) END IF IF (aa <= xtol*ABS(cs)) EXIT ak = ak + 1.0_wp END DO IF (nd == 1) cs = cs + (-LOG(z) + euler) IF (kode == 2) cs = cs * EXP(z) cy(1) = cs ct = cs emz = cone IF (m /= 1) THEN cy(ind) = cs ak = ks IF (kode == 1) emz = EXP(-z) SELECT CASE ( icase ) CASE ( 1) GO TO 140 CASE ( 2) GO TO 160 END SELECT END IF IF (icase == 2) RETURN IF (kode == 1) emz = EXP(-z) cy(1) = (emz-cs) / z RETURN END IF !----------------------------------------------------------------------- ! BACKWARD RECURSIVE MILLER ALGORITHM FOR ! E(N,Z)=EXP(-Z)*(Z**(N-1))*U(N,N,Z) ! WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO ABS(Z) ! U(A,B,Z) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION !----------------------------------------------------------------------- emz = cone IF (kode /= 2) THEN !----------------------------------------------------------------------- ! SCALE NEAR EXPONENT EXTREMES ON KODE=1 !----------------------------------------------------------------------- IF (x >= 0.0_wp) THEN at = x + (n+m-1) ct = CMPLX(at, y, KIND=wp) aa = ABS(ct) at = x + LOG(aa) IF (at > elim) GO TO 30 kflag = 1 ELSE at = x IF (at < (-elim)) GO TO 180 kflag = 2 END IF IF (ABS(at) >= alim) THEN tola = EXP(alim-elim) rtola = 1.0_wp / tola IF (kflag /= 2) THEN scle = rtola rscle = tola ELSE scle = tola rscle = rtola END IF emz = scle END IF emz = emz * EXP(-z) END IF iz = az + 0.5_wp kn = n + m - 1 IF (kn <= iz) GO TO 50 IF (n < iz .AND. iz < kn) GO TO 80 IF (n >= iz) GO TO 70 ierr = 7 RETURN 30 ierr = 2 cy(1:m) = czero RETURN 50 icase = 1 ks = kn ml = m - 1 mu = -1 ind = m IF (kn > 1) GO TO 90 60 ks = 2 icase = 3 GO TO 90 70 icase = 2 ind = 1 ks = n mu = m - 1 IF (n > 1) GO TO 90 IF (kn == 1) GO TO 60 iz = 2 80 icase = 1 ks = iz ml = iz - n ind = ml + 1 mu = kn - iz 90 ik = ks / 2 ah = ik jset = 1 + ks - 2 * ik !----------------------------------------------------------------------- ! START COMPUTATION FOR ! CYY(1) = C*U( A , A ,Z) JSET=1 ! CYY(1) = C*U(A+1,A+1,Z) JSET=2 ! FOR AN EVEN INTEGER A. !----------------------------------------------------------------------- ic = 0 aa = ah + ah caa = aa aam = aa - 1.0_wp aams = aam * aam tz = z + z fz = tz + tz ak = ah xtol = tol ct = aams + fz * ah cak = z + caa aem = (ak+1.0_wp) / xtol aem = aem / ABS(cak) bk = aa ck = ah * ah !----------------------------------------------------------------------- ! FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD ! RECURSION !----------------------------------------------------------------------- cp1 = czero cp2 = cone 100 ic = ic + 1 IF (ic <= icdim) THEN ak = ak + 1.0_wp ck = ck + 1.0_wp at = bk / (bk+ak+ck) bk = bk + ak + ak cat = at ca(ic) = at bt = 1.0_wp / (ak + 1.0_wp) cbt = (ak + ak + z) * bt cb(ic) = cbt cpt = cp2 cp2 = cbt * cp2 - cat * cp1 cp1 = cpt ct = ct + fz aem = aem * at bt = ABS(ct) ERR = aem / SQRT(bt) ap1 = ABS(cp1) IF (ERR*(ak+1.0_wp)/ap1 > ap1) GO TO 100 ELSE !----------------------------------------------------------------------- ! CONTINUE FORWARD RECURRENCE UNINDEXED WHEN IC EXCEEDS ICDIM !----------------------------------------------------------------------- ic = ic - 1 110 ic = ic + 1 IF (ic > icmax) GO TO 190 ak = ak + 1.0_wp ck = ck + 1.0_wp at = bk / (bk+ak+ck) bk = bk + ak + ak cat = at bt = 1.0_wp / (ak+1.0_wp) cbt = (ak+ak + z) * bt cpt = cp2 cp2 = cbt * cp2 - cat * cp1 cp1 = cpt ct = ct + fz aem = aem * at bt = ABS(ct) ERR = aem / SQRT(bt) ap1 = ABS(cp1) IF (ERR*(ak+1.0_wp)/ap1 > ap1) GO TO 110 END IF fc = ic at = ((fc+1.0_wp)/(ak+1.0_wp)) * ((ak+ah)/(ak+1.0_wp)) cat = at * SQRT(ct/(ct+fz)) cy2 = cat * (cp1/cp2) cy1 = cone !----------------------------------------------------------------------- ! BACKWARD RECURRENCE FOR ! CY1=C*U( A ,A,Z) ! CY2=C*(A/(1+A/2))*U(A+1,A,Z) !----------------------------------------------------------------------- IF (ic > icdim) THEN !----------------------------------------------------------------------- ! BACKWARD RECUR UNINDEXED WHEN IC EXCEEDS ICDIM !----------------------------------------------------------------------- bt = aa + x 120 ak = ah + fc bk = ak + 1.0_wp ck = aam + fc dk = bt + fc + fc at = (ak/fc) * (bk/ck) cbt = CMPLX(dk/bk, y/bk, KIND=wp) cpt = cy1 cy1 = (cbt*cy1 - cy2) * at cy2 = cpt fc = fc - 1.0_wp ic = ic - 1 IF (ic > icdim) GO TO 120 END IF ict = ic DO k = 1, ict at = 1.0_wp / ca(ic) cpt = cy1 cy1 = (cb(ic)*cy1 - cy2) * at cy2 = cpt ic = ic - 1 END DO !----------------------------------------------------------------------- ! THE CONTIGUOUS RELATION ! Z*U(B,C+1,Z)=(C-B)*U(B,C,Z)+U(B-1,C,Z) ! WITH B=A+1 , C=A IS USED FOR ! CYY(2) = C * U(A+1,A+1,Z) ! Z IS INCORPORATED INTO THE NORMALIZING RELATION FOR CNORM. !----------------------------------------------------------------------- cpt = cy2 / cy1 cnorm = cone - cpt * (ah + 1.0_wp) / aa cyy(1) = cone / (cnorm*caa+z) cyy(2) = cnorm * cyy(1) IF (icase /= 3) THEN ct = emz * cyy(jset) cy(ind) = ct * rscle cs = ct IF (m == 1) RETURN ak = ks SELECT CASE ( icase ) CASE ( 1) GO TO 140 CASE ( 2) GO TO 160 END SELECT END IF !----------------------------------------------------------------------- ! RECURSION SECTION N*E(N+1,Z) + Z*E(N,Z)=EMZ !----------------------------------------------------------------------- ct = emz * (cone - cyy(1)) / z cy(1) = ct * rscle RETURN 140 caa = ak tz = cone / z k = ind - 1 DO i = 1, ml caa = caa - cone ct = (emz-caa*ct) * tz cy(k) = ct * rscle k = k - 1 END DO IF (mu <= 0) RETURN ak = ks 160 k = ind DO i = 1, mu cs = (emz - z*cs) / ak cy(k+1) = cs * rscle ak = ak + 1.0_wp k = k + 1 END DO RETURN 180 ierr = 3 RETURN 190 ierr = 6 RETURN END SUBROUTINE cexenz SUBROUTINE cacexi(z, nu, kode, n, y, ierr, yb, rbry, tol, elim, alim, icdim, ca, cb) ! REFER TO CEXINT ! ! CACEXI COMPUTES THE ANALYTIC CONTINUATION OF THE EXPONENTIAL INTEGRAL ! FOR X < 0 AND -YB < Y < YB BY TAYLOR SERIES IN INCREMENTS OF HALF A UNIT. ! THE CONTINUATION PROCEEDS VERTICALLY DOWN FROM ZB=CMPLX(X,YB) TO ! Z=CMPLX(X,Y) FOR 0.0 <= Y < YB. CONJUGATION IS USED FOR Y < 0.0E0. ! ! ROUTINES CALLED CEXENZ COMPLEX(wp), INTENT(IN) :: z INTEGER, INTENT(IN) :: nu INTEGER, INTENT(IN) :: kode INTEGER, INTENT(IN) :: n COMPLEX(wp), INTENT(OUT) :: y(n) INTEGER, INTENT(OUT) :: ierr REAL(wp), INTENT(IN OUT) :: yb REAL(wp), INTENT(IN OUT) :: rbry REAL(wp), INTENT(IN) :: tol REAL(wp), INTENT(IN) :: elim REAL(wp), INTENT(IN) :: alim INTEGER, INTENT(IN) :: icdim REAL(wp), INTENT(IN OUT) :: ca(icdim) COMPLEX(wp), INTENT(IN OUT) :: cb(icdim) INTEGER :: i, iaz, il, is, iy, k, kb, kl, kmax, ks, kyb, nb, nflg, nub REAL(wp) :: az, del, fk, rtola, tola, yt, zi, zid, zr, & xtol, rzi, atrm, fj, rw, rq(64), asum, htol COMPLEX(wp) :: yy(1), cex, sum, trm, zz, zp(64), zt, scle, rscle, trms, zw, cezt COMPLEX(wp), PARAMETER :: cone = (1.0_wp, 0.0_wp) scle = cone rscle = cone zr = REAL(z, KIND=wp) zi = AIMAG(z) zid = zi kyb = 0 IF (zi < 0.0_wp) zid = -zid az = ABS(z) iaz = az + 0.5_wp !----------------------------------------------------------------------- ! SET ORDER NUB=MIN(N+M-1,INT(ABS(Z)+0.5)), GENERATE REMAINDER ! OF THE SEQUENCE AFTER THE COMPUTATION OF E(NUB,Z) !----------------------------------------------------------------------- IF (nu < iaz) THEN IF (nu+n-1 <= iaz) GO TO 10 nub = iaz nb = 0 nflg = 3 kb = nub - nu + 1 ks = kb + 1 kl = n is = 1 il = kb - 1 GO TO 30 END IF nub = MAX(iaz, 1) nb = nu - nub nflg = 1 kb = 1 ks = 2 kl = n GO TO 30 10 nub = nu + n - 1 nb = 0 nflg = 2 is = 2 il = n kb = n GO TO 30 !----------------------------------------------------------------------- ! SET PARAMETERS FOR ANALYTIC CONTINUATION FROM Y=YB INTO THE REGION ! 0 <= ZID < YB. !----------------------------------------------------------------------- 20 yb = yb + 0.5_wp kyb = kyb + 1 IF (kyb > 10) RETURN ierr = 0 30 del = yb - zid !----------------------------------------------------------------------- ! MAKE DEL LARGE ENOUGH TO AVOID UNDERFLOW IN GENERATION OF POWERS !----------------------------------------------------------------------- IF (ABS(del) <= 1.0E-4_wp) THEN yb = yb + 1.0E-4_wp del = yb - zid END IF htol = 0.125_wp * tol zz = CMPLX(zr, yb, KIND=wp) CALL cexenz(zz, nub, 2, 1, yy, ierr, rbry, htol, elim, alim, icdim, ca, cb) IF (ierr == 6) GO TO 20 !----------------------------------------------------------------------- ! ANALYTIC CONTINUATION VIA TAYLOR SERIES FOR ORDER NUB !----------------------------------------------------------------------- iy = del + del yt = del / (iy+1) sum = yy(1) htol = 0.25_wp * tol trm = sum cezt = CMPLX(COS(yt), -SIN(yt), KIND=wp) zw = cone zp(1) = cone fk = 1.0_wp fj = nub - 1 zt = cone / zz !----------------------------------------------------------------------- ! TERMS OF THE SERIES TRM=E(NUB-K,ZZ)*(YT**K)/K!, K=0,1,... ARE ! GENERATED BY BACKWARD RECURRENCE. E IS SCALED BY CEXP(ZZ). !----------------------------------------------------------------------- DO k = 2, 64 rw = yt / fk zw = zw * CMPLX(0.0_wp, rw, KIND=wp) rw = fj * rw trm = (zw - CMPLX(0.0_wp, rw, KIND=wp)*trm) * zt sum = sum + trm zp(k) = zw rq(k) = rw asum = ABS(sum) atrm = ABS(trm) IF (atrm < htol*asum) GO TO 50 fk = fk + 1.0_wp fj = fj - 1.0_wp END DO k = 64 50 kmax = k sum = sum * cezt IF (iy /= 0) THEN DO i = 1, iy rzi = (iy-i+1) * yt + zid zz = CMPLX(zr, rzi, KIND=wp) zt = cone / zz trm = sum DO k = 2, kmax trm = (zp(k) - CMPLX(0.0_wp, rq(k), KIND=wp)*trm) * zt sum = sum + trm END DO atrm = ABS(trm) asum = ABS(sum) xtol = htol * asum IF (atrm >= xtol) THEN IF (kmax < 64) THEN kmax = kmax + 1 DO k = kmax, 64 rw = yt / fk zw = zw * CMPLX(0.0_wp, rw, KIND=wp) rw = fj * rw trm = (zw - CMPLX(0.0_wp, rw, KIND=wp)*trm) * zt sum = sum + trm zp(k) = zw rq(k) = rw atrm = ABS(trm) IF (atrm < xtol) GO TO 80 fk = fk + 1.0_wp fj = fj - 1.0_wp END DO k = 64 80 kmax = k END IF END IF sum = sum * cezt END DO END IF !----------------------------------------------------------------------- ! FORM SEQUENCE UP OR DOWN FROM ORDER NUB !----------------------------------------------------------------------- IF (zi < 0.0_wp) sum = CONJG(sum) cex = cone !----------------------------------------------------------------------- ! SCALE NEAR OVERFLOW LIMIT ON KODE=1 !----------------------------------------------------------------------- IF (kode /= 2) THEN IF (ABS(zr) >= alim) THEN IF (ABS(zr) > elim) GO TO 130 tola = EXP(alim-elim) rtola = 1.0_wp / tola scle = tola rscle = rtola END IF cex = scle * EXP(-z) END IF trm = sum * cex y(kb) = trm * rscle trms = trm IF (n == 1 .AND. nflg /= 1) RETURN IF (nflg /= 2) THEN fk = nub IF (nflg == 1 .AND. nb /= 0) THEN DO k = 1, nb trm = (cex - z*trm) / fk fk = fk + 1.0_wp END DO END IF y(kb) = trm * rscle trms = trm IF (n == 1) RETURN DO k = ks, kl trm = (cex - z*trm) / fk y(k) = trm * rscle fk = fk + 1.0_wp END DO IF (nflg == 1) RETURN END IF k = kb - 1 fk = nub - 1 zt = cone / z DO i = is, il trms = (cex - fk*trms) * zt y(k) = trms * rscle k = k - 1 fk = fk - 1.0_wp END DO RETURN 130 ierr = 3 RETURN END SUBROUTINE cacexi FUNCTION psixn(n) RESULT(fn_val) !! CALGO 683 Digamma function \(\psi_0(n)\). ! !! \(\lbrace n \in \mathbb{Z} \mid n \geq 1 \rbrace\) ! ! REFER TO CEXENZ ! ! THIS SUBROUTINE RETURNS VALUES OF PSI(X)=DERIVATIVE OF LOG GAMMA(X), ! X > 0.0 AT INTEGER ARGUMENTS. A TABLE LOOK-UP IS PERFORMED FOR N <= 100, ! AND THE ASYMPTOTIC EXPANSION IS EVALUATED FOR N > 100. INTEGER, INTENT(IN) :: n REAL(wp) :: fn_val INTEGER :: k REAL(wp) :: ax, fn, rfn2, trm, s, wdtol !----------------------------------------------------------------------- ! PSIXN(N), N = 1,100 !----------------------------------------------------------------------- REAL(wp), PARAMETER :: c(100) = (/ -5.77215664901532861e-1_wp, & 4.22784335098467139e-1_wp, 9.22784335098467139e-1_wp, 1.25611766843180047_wp, & 1.50611766843180047_wp, 1.70611766843180047_wp, 1.87278433509846714_wp, & 2.01564147795561000_wp, 2.14064147795561000_wp, 2.25175258906672111_wp, & 2.35175258906672111_wp, 2.44266167997581202_wp, 2.52599501330914535_wp, & 2.60291809023222227_wp, 2.67434666166079370_wp, 2.74101332832746037_wp, & 2.80351332832746037_wp, 2.86233685773922507_wp, 2.91789241329478063_wp, & 2.97052399224214905_wp, 3.02052399224214905_wp, 3.06814303986119667_wp, & 3.11359758531574212_wp, 3.15707584618530734_wp, 3.19874251285197401_wp, & 3.23874251285197401_wp, 3.27720405131351247_wp, 3.31424108835054951_wp, & 3.34995537406483522_wp, 3.38443813268552488_wp, 3.41777146601885821_wp, & 3.45002953053498724_wp, 3.48127953053498724_wp, 3.51158256083801755_wp, & 3.54099432554389990_wp, 3.56956575411532847_wp, 3.59734353189310625_wp, & 3.62437055892013327_wp, 3.65068634839381748_wp, 3.67632737403484313_wp, & 3.70132737403484313_wp, 3.72571761793728215_wp, 3.74952714174680596_wp, & 3.77278295570029433_wp, 3.79551022842756706_wp, 3.81773245064978928_wp, & 3.83947158108457189_wp, 3.86074817682925274_wp, 3.88158151016258607_wp, & 3.90198967342789220_wp, 3.92198967342789220_wp, 3.94159751656514710_wp, & 3.96082828579591633_wp, 3.97969621032421822_wp, 3.99821472884273674_wp, & 4.01639654702455492_wp, 4.03425368988169777_wp, 4.05179754953082058_wp, & 4.06903892884116541_wp, 4.08598808138353829_wp, 4.10265474805020496_wp, & 4.11904819067315578_wp, 4.13517722293122029_wp, 4.15105023880423617_wp, & 4.16667523880423617_wp, 4.18205985418885155_wp, 4.19721136934036670_wp, & 4.21213674247469506_wp, 4.22684262482763624_wp, 4.24133537845082464_wp, & 4.25562109273653893_wp, 4.26970559977879245_wp, 4.28359448866768134_wp, & 4.29729311880466764_wp, 4.31080663231818115_wp, 4.32413996565151449_wp, & 4.33729786038835659_wp, 4.35028487337536958_wp, 4.36310538619588240_wp, & 4.37576361404398366_wp, 4.38826361404398366_wp, 4.40060929305632934_wp, & 4.41280441500754886_wp, 4.42485260777863319_wp, 4.43675736968339510_wp, & 4.44852207556574804_wp, 4.46014998254249223_wp, 4.47164423541605544_wp, & 4.48300787177969181_wp, 4.49424382683587158_wp, 4.50535493794698269_wp, & 4.51634394893599368_wp, 4.52721351415338499_wp, 4.53796620232542800_wp, & 4.54860450019776842_wp, 4.55913081598724211_wp, 4.56954748265390877_wp, & 4.57985676100442424_wp, 4.59006084263707730_wp, 4.60016185273808740_wp /) !----------------------------------------------------------------------- ! COEFFICIENTS OF ASYMPTOTIC EXPANSION !----------------------------------------------------------------------- REAL(wp), PARAMETER :: b(6) = (/ & 8.33333333333333333e-2_wp, -8.33333333333333333e-3_wp, & 3.96825396825396825e-3_wp, -4.16666666666666666e-3_wp, & 7.57575757575757576e-3_wp, -2.10927960927960928e-2_wp & /) IF (n <= 100) THEN fn_val = c(n) RETURN END IF wdtol = MAX(EPSILON(0.0_wp), 1.0e-18_wp) fn = n ax = 1.0_wp s = -0.5_wp / fn IF (ABS(s) > wdtol) THEN rfn2 = 1.0_wp / (fn*fn) DO k = 1, 6 ax = ax * rfn2 trm = -b(k) * ax IF (ABS(trm) < wdtol) EXIT s = s + trm END DO END IF fn_val = s + LOG(fn) RETURN END FUNCTION psixn END MODULE calgo_683