! Licensed under the Computer Physics Communications (CPC) User License ! Copyright © 2018 N. Michel ! See ColSpecF LICENSE file for details. module cpc_michel !* # CPC Michel ! CPC Algorithm AEAE_v1_0. ! ! Procedures: ! ! - `hyp_2f1`: Gauss hypergeometric function \({}_2F_1(a, b; c; z)\) ! ! ## Author ! N. Michel ! ! ## History ! - 2007-12-17 - N. Michel ! - Original code. ! - 2008-04-01 - N. Michel ! - F90 code distributed by CPC: ! ! - 2025-07-12 - Rodrigo Castro (GitHub: rodpcastro) ! - Placed all procedures in the module `cpc_michel` and removed references to the ! old module `HYP_2F1_MODULE`. ! - Removed unnecessary declarations of function type. ! - Positioned original docstrings inside their respective procedures. ! - Replaced `PR` (real and complex precision) by `wp` (CSF working precision). ! - Replaced `IPR` (integer precision) by `i4` (CSF 4-byte integer). ! - 2025-07-14 - Rodrigo Castro (GitHub: rodpcastro) ! - Replaced `EPS15` by `eps_wp` (CSF working precision machine epsilon). ! ! ## References ! 1. N. Michel and M. V. Stoitsov. 2008. Fast computation of the Gauss hypergeometric ! function with all its parameters complex with application to the Pöschl-Teller- ! Ginocchio potential wave functions. Computer Physics Communications 178, 7 (April !* 2008), 535–551. use csf_kinds, only: wp, i4 use csf_numerror, only: eps_wp implicit none private public :: hyp_2f1 REAL(wp) :: ZERO=0.0_wp,ONE=1.0_wp,TWO=2.0_wp,HALF=0.50_wp REAL(wp) :: M_PI=3.14159265358979323846_wp REAL(wp) :: M_PI_2=1.57079632679489661923_wp REAL(wp) :: M_1_PI=0.31830988618379067154_wp CONTAINS FUNCTION INF_NORM(Z) COMPLEX(wp),INTENT(IN) :: Z REAL(wp) :: INF_NORM INF_NORM=MAX(ABS(REAL(Z,wp)),ABS(AIMAG(Z))) RETURN END FUNCTION INF_NORM FUNCTION TANZ(Z) COMPLEX(wp),INTENT(IN) :: Z COMPLEX(wp) :: TANZ TANZ=SIN(Z)/COS(Z) RETURN END FUNCTION TANZ FUNCTION LOG1P(Z) ! log1p(z) = log(1.0 + z) ! ! Reference ! --------- ! 1. N. J. Higham. 1996. Accuracy and Stability ! of Numerical Algorithms. SIAM, Philadelphia. COMPLEX(wp),INTENT(IN) :: Z REAL(wp) :: X,XP1,LOG1P_X REAL(wp) :: Y,YX,YX2,YX2P1,LOG1P_YX2 REAL(wp) :: RE_LOG1P,IM_LOG1P COMPLEX(wp) :: LOG1P IF(INF_NORM(Z).LT.ONE) THEN X = REAL(Z,wp); XP1 = X+ONE IF(XP1.EQ.ONE) THEN LOG1P_X = X ELSE LOG1P_X = LOG(XP1)*X/(XP1-ONE) ENDIF Y = AIMAG(Z) YX = Y/XP1; YX2 = YX*YX; YX2P1 = YX2+ONE IF(YX2P1.EQ.ONE) THEN LOG1P_YX2 = YX2 ELSE LOG1P_YX2 = LOG(YX2P1)*YX2/(YX2P1-ONE) ENDIF RE_LOG1P = LOG1P_X + HALF*LOG1P_YX2 IM_LOG1P = ATAN2(Y,XP1) LOG1P = CMPLX(RE_LOG1P,IM_LOG1P,wp) RETURN ELSE LOG1P=LOG(ONE+Z) RETURN ENDIF END FUNCTION LOG1P FUNCTION EXPM1(Z) ! expm1(z) = exp(z) - 1.0 ! ! Reference ! --------- ! 1. N. J. Higham. 1996. Accuracy and Stability ! of Numerical Algorithms. SIAM, Philadelphia. COMPLEX(wp),INTENT(IN) :: Z REAL(wp) :: X,EXPM1_X,EXP_X,Y,SIN_HALF_Y REAL(wp) :: RE_EXPM1,IM_EXPM1 COMPLEX(wp) :: EXPM1 IF(INF_NORM(Z).LT.ONE) THEN X = REAL(Z,wp); EXP_X = EXP(X) Y = AIMAG(Z); SIN_HALF_Y=SIN(HALF*Y) IF(EXP_X.EQ.ONE) THEN EXPM1_X = X ELSE EXPM1_X = (EXP_X-ONE)*X/LOG(EXP_X) ENDIF RE_EXPM1 = EXPM1_X-TWO*EXP_X*SIN_HALF_Y*SIN_HALF_Y IM_EXPM1 = EXP_X*SIN(Y) EXPM1 = CMPLX(RE_EXPM1,IM_EXPM1,wp) RETURN ELSE EXPM1=EXP(Z)-ONE RETURN ENDIF END FUNCTION EXPM1 RECURSIVE FUNCTION LOG_GAMMA(Z) RESULT(RES) ! Logarithm of Gamma[z] and Gamma inverse function ! ! For log[Gamma[z]],if z is not finite ! or is a negative integer, the program ! returns an error message and stops. ! The Lanczos method is used. Precision : ~ 1E-15 ! The method works for Re[z]>0.5 . ! If Re[z]<=0.5, one uses the formula Gamma[z].Gamma[1-z]=Pi/sin(Pi.z) ! log[sin(Pi.z)] is calculated with the Kolbig method ! (K.S. Kolbig, Comp. Phys. Comm., Vol. 4, p.221(1972)): ! If z=x+iy and y>=0, log[sin(Pi.z)]=log[sin(Pi.eps)]-i.Pi.n, ! with z=n+eps so 0<=Re[eps]< 1 and n integer. ! If y>110, log[sin(Pi.z)]=-i.Pi.z+log[0.5]+i.Pi/2 ! numerically so that no overflow can occur. ! If z=x+iy and y< 0, log[Gamma(z)]=[log[Gamma(z*)]]*, ! so that one can use the previous formula with z*. ! ! For Gamma inverse, Lanczos method is also used ! with Euler reflection formula. ! sin (Pi.z) is calculated as sin (Pi.(z-n)) ! to avoid inaccuracy with z = n + eps ! with n integer and |eps| as small as possible. ! ! ! Variables: ! ---------- ! x,y: Re[z], Im[z] ! log_sqrt_2Pi,log_Pi : log[sqrt(2.Pi)], log(Pi). ! sum : Rational function in the Lanczos method ! log_Gamma_z : log[Gamma(z)] value. ! c : table containing the fifteen coefficients in the expansion ! used in the Lanczos method. ! eps,n : z=n+eps so 0<=Re[eps]< 1 and n integer for Log[Gamma]. ! z=n+eps and n integer ! so |eps| is as small as possible for Gamma_inv. ! log_const : log[0.5]+i.Pi/2 ! g : coefficient used in the Lanczos formula. It is here 607/128. ! z,z_m_0p5,z_p_g_m0p5,zm1 : argument of the Gamma function, ! z-0.5, z-0.5+g, z-1 ! res: returned value !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: Z INTEGER(i4) :: N,I REAL(wp) :: X,Y,LOG_SQRT_2PI,G,LOG_PI,M_LN2,C(0:14) COMPLEX(wp) :: GAMMA_SUM,Z_M_0P5,Z_P_G_M0P5,ZM1 COMPLEX(wp) :: LOG_CONST,I_PI,EPS,LOG_SIN_PI_Z,RES ! M_LN2=0.69314718055994530942_wp; X=REAL(Z,wp); Y=AIMAG(Z) IF((Z.EQ.NINT(X)).AND.(X.LE.ZERO)) & STOP 'Z IS NEGATIVE INTEGER IN LOG_GAMMA' IF(X.GE.HALF) THEN LOG_SQRT_2PI=0.91893853320467274177_wp; G=4.7421875_wp Z_M_0P5=Z-HALF; Z_P_G_M0P5=Z_M_0P5+G; ZM1=Z-ONE C=(/ 0.99999999999999709182_wp,57.156235665862923517_wp, & -59.597960355475491248_wp, 14.136097974741747174_wp, & -0.49191381609762019978_wp, 0.33994649984811888699e-4_wp, & 0.46523628927048575665e-4_wp, -0.98374475304879564677e-4_wp, & 0.15808870322491248884e-3_wp, -0.21026444172410488319e-3_wp, & 0.21743961811521264320e-3_wp, -0.16431810653676389022e-3_wp, & 0.84418223983852743293e-4_wp, -0.26190838401581408670e-4_wp, & 0.36899182659531622704e-5_wp /) GAMMA_SUM=C(0) DO I=1,14 GAMMA_SUM=GAMMA_SUM+C(I)/(ZM1+I) ENDDO RES=LOG_SQRT_2PI+LOG(GAMMA_SUM)+Z_M_0P5*LOG(Z_P_G_M0P5) & -Z_P_G_M0P5 RETURN ELSE IF(Y.GE.ZERO) THEN IF(X.LT.NINT(X)) THEN N=NINT(X)-1 ELSE N=NINT(X) ENDIF LOG_PI=1.1447298858494002_wp LOG_CONST=CMPLX(-M_LN2,M_PI_2,wp); I_PI=CMPLX(ZERO,M_PI,wp) EPS=Z-N IF(Y.GT.110.0_wp) THEN LOG_SIN_PI_Z=-I_PI*Z+LOG_CONST ELSE LOG_SIN_PI_Z=LOG(SIN(M_PI*EPS))-I_PI*N ENDIF RES=LOG_PI-LOG_SIN_PI_Z-LOG_GAMMA(ONE-Z); RETURN ELSE RES=CONJG(LOG_GAMMA(CONJG(Z))) RETURN ENDIF END FUNCTION LOG_GAMMA RECURSIVE FUNCTION GAMMA_INV(Z) RESULT(RES) ! Inverse of the Gamma function [1/Gamma](z) ! ! It is calculated with the Lanczos method for Re[z] >= 0.5 ! and is precise up to 10^{-15}. ! If Re[z] <= 0.5, one uses the formula ! Gamma[z].Gamma[1-z] = Pi/sin (Pi.z). ! sin (Pi.z) is calculated as sin (Pi.(z-n)) to avoid inaccuracy, ! with z = n + eps with n integer and |eps| as small as possible. ! ! Variables ! --------- ! z : argument of the function ! x: Re[z] ! eps,n : z = n + eps with n integer and |eps| as small as possible. ! res: returned value !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: Z INTEGER(i4) :: N,I REAL(wp) :: X,LOG_SQRT_2PI,G,C(0:14) COMPLEX(wp) :: RES,GAMMA_SUM,Z_M_0P5,Z_P_G_M0P5,ZM1,EPS ! X=REAL(Z,wp) IF(X.GE.HALF) THEN LOG_SQRT_2PI=0.91893853320467274177_wp; G=4.7421875_wp Z_M_0P5=Z-HALF; Z_P_G_M0P5=Z_M_0P5+G; ZM1=Z-ONE C=(/ 0.99999999999999709182_wp,57.156235665862923517_wp, & -59.597960355475491248_wp, 14.136097974741747174_wp, & -0.49191381609762019978_wp, 0.33994649984811888699e-4_wp, & 0.46523628927048575665e-4_wp, -0.98374475304879564677e-4_wp, & 0.15808870322491248884e-3_wp, -0.21026444172410488319e-3_wp, & 0.21743961811521264320e-3_wp, -0.16431810653676389022e-3_wp, & 0.84418223983852743293e-4_wp, -0.26190838401581408670e-4_wp, & 0.36899182659531622704e-5_wp /) GAMMA_SUM=C(0) DO I=1,14 GAMMA_SUM=GAMMA_SUM+C(I)/(ZM1+I); ENDDO RES=EXP(Z_P_G_M0P5-Z_M_0P5*LOG(Z_P_G_M0P5)-LOG_SQRT_2PI) & /GAMMA_SUM RETURN ELSE X=REAL(Z,wp); N=NINT(X) EPS=Z-N IF(MOD(N,2).EQ.0) THEN RES=SIN(M_PI*EPS)*M_1_PI/GAMMA_INV (ONE-Z) RETURN ELSE RES=-SIN(M_PI*EPS)*M_1_PI/GAMMA_INV (ONE-Z) RETURN ENDIF ENDIF END FUNCTION GAMMA_INV RECURSIVE FUNCTION GAMMA_RATIO_DIFF_SMALL_EPS(Z,EPS) RESULT(RES) ! Calculation of H(z,eps) = [Gamma(z+eps)/Gamma(z) - 1]/eps, with e and ! z complex so z,z+eps are not negative integers and 0 <= |eps|oo < 0.1 ! ! The function H(z,eps) = [Gamma(z+eps)/Gamma(z) - 1]/e is calculated ! here with the Lanczos method. ! For the Lanczos method, the gamma parameter, denoted as g, ! is 4.7421875 and one uses a sum of 15 numbers with the table c[15], ! so that it is precise up to machine accuracy. ! The H(z,eps) function is used in formulas occuring in1-z and 1/z ! transformations (see Comp. Phys. Comm. paper). ! ! One must have z and z+eps not negative integers as otherwise ! it is clearly not defined. ! As this function is meant to be precise for small |eps|oo, ! one has to have 0 <= |eps|oo < 0.1 . ! Indeed, a direct implementation of H(z,eps) with Gamma_inv or ! log_Gamma for |eps|oo >= 0.1 is numerically stable. ! The returned function has full numerical accuracy ! even if |eps|oo is very small. ! ! eps not equal to zero ! --------------------- ! If Re(z) >= 0.5 or Re(z+eps) >= 0.5, one clearly has Re(z) > 0.4 ! and Re(z+eps) > 0.4, ! so that the Lanczos summation can be used for both Gamma(z) ! and Gamma(z+eps). ! One then has: ! log[Gamma(z+eps)/Gamma(z)] = ! (z-0.5) log1p[eps/(z+g-0.5)] + eps log(z+g-0.5+eps) - eps ! + log1p[-eps \sum_{i=1}^{14} c[i]/((z-1+i)(z-1+i+eps)) ! / (c[0] + \sum_{i=1}^{14} c[i]/(z-1+i))] ! H(z,eps) = expm1[log[Gamma(z+eps)/Gamma(z)]]/eps . ! ! If Re(z) < 0.5 and Re(z+eps) < 0.5, ! Euler reflection formula is used for both Gamma(z) and Gamma(z+eps). ! One then has: ! H(z+eps,-eps) = [cos(pi.eps) + sin(pi.eps)/tan(pi(z-n))].H(1-z,-eps) ! + (2/eps).sin^2(eps.pi/2) - sin(pi.eps)/(eps.tan(pi.(z-n))) ! H(1-z,-eps) is calculated with the Lanczos summation ! as Re(1-z) >= 0.5 and Re(1-z-eps) >= 0.5 . ! z-n is used in tan(pi.z) instead of z to avoid inaccuracies ! due the finite number of digits of pi. ! H(z,eps) = H(z+eps,-eps)/(1 - eps.H(z+eps,-eps)) ! provides the final result. ! ! eps equal to zero ! ----------------- ! It is obtained with the previous case and eps -> 0 : ! If Re(z) >= 0.5, one has: ! H(z,eps) = (z-0.5)/(z+g-0.5) + log(z+g-0.5) - 1 - ! \sum_{i=1}^{14} c[i]/((z-1+i)^2)/(c[0]+\sum_{i=1}^{14} c[i]/(z-1+i)) ! ! If Re(z) < 0.5, one has: ! H(z,0) = H(1-z,0) - pi/tan(pi.(z-n)) ! ! Variables ! --------- ! z,eps: input variables of the function H(z,eps) ! g,c[15]: double and table of 15 doubles defining the Lanczos sum ! so that it provides the Gamma function ! precise up to machine accuracy. ! eps_pz,z_m_0p5,z_pg_m0p5,eps_pz_pg_m0p5,zm1,zm1_p_eps: ! z+eps,z-0.5,z+g-0.5,z+eps+g-0.5,z-1,z-1+eps ! x,eps_px: real parts of z and z+eps. ! n,m: closest integer ot the real part of z, same for z+eps. ! sum_num,sum_den: \sum_{i=1}^{14} c[i]/((z-1+i)(z-1+i+eps)) ! and (c[0] + \sum_{i=1}^{14} c[i]/(z-1+i)). ! They appear respectively as numerator and denominator in formulas. ! Pi_eps,term,T1_eps_z: pi.eps, sin (pi.eps)/tan(pi.(z-n)), ! [cos(pi.eps) + sin(pi.eps)/tan(pi(z-n))].H(1-z,-eps) ! sin_Pi_2_eps,T2_eps_z,T_eps_z: sin^2(eps.pi/2), ! (2/eps).sin^2(eps.pi/2) - sin(pi.eps)/(eps.tan(pi.(z-n))), ! H(z+eps,-eps) ! res: returned value !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: Z,EPS INTEGER(i4) :: N,M,I REAL(wp) :: G,X,EPS_PX,C(0:14) COMPLEX(wp) :: RES,SUM_NUM,SUM_DEN COMPLEX(wp) :: EPS_PZ,Z_M_0P5,Z_PG_M0P5,EPS_PZ_PG_M0P5,ZM1 COMPLEX(wp) :: CI_ZM1_PI_INV,PI_EPS,TT,T1_EPS_Z,SIN_PI_2_EPS COMPLEX(wp) :: ZM1_P_EPS,T2_EPS_Z,T_EPS_Z ! G=4.74218750_wp IF(INF_NORM(EPS).GT.0.1_wp) & STOP 'ONE MUST HAVE |EPS|< 0.1 IN GAMMA_RATIO_DIFF_SMALL_EPS' EPS_PZ=Z+EPS; Z_M_0P5=Z-HALF; Z_PG_M0P5=Z_M_0P5+G EPS_PZ_PG_M0P5=Z_PG_M0P5+EPS; ZM1=Z-ONE; ZM1_P_EPS=ZM1+EPS X=REAL(Z,wp); EPS_PX=REAL(EPS_PZ,wp); N=NINT(X); M=NINT(EPS_PX) IF((Z.EQ.N).AND.(N.LE.0)) THEN STOP 'Z IS NEGATIVE INTEGER IN GAMMA_RATIO_DIFF_SMALL_EPS' ENDIF IF((EPS_PZ.EQ.M).AND.(M.LE.0)) THEN STOP 'Z+EPS IS NEGATIVE INTEGER IN GAMMA_RATIO_DIFF_SMALL_EPS' ENDIF C=(/ 0.99999999999999709182_wp,57.156235665862923517_wp, & -59.597960355475491248_wp,14.136097974741747174_wp, & -0.49191381609762019978_wp,0.33994649984811888699e-4_wp, & 0.46523628927048575665e-4_wp,-0.98374475304879564677e-4_wp, & 0.15808870322491248884e-3_wp,-0.21026444172410488319e-3_wp, & 0.21743961811521264320e-3_wp,-0.16431810653676389022e-3_wp, & 0.84418223983852743293e-4_wp,-0.26190838401581408670e-4_wp, & 0.36899182659531622704e-5_wp /) IF((X.GE.HALF).OR.(EPS_PX.GE.HALF)) THEN SUM_NUM=ZERO;SUM_DEN=C(0) DO I=1,14 CI_ZM1_PI_INV=C(I)/(ZM1+I) SUM_NUM=SUM_NUM+CI_ZM1_PI_INV/(ZM1_P_EPS+I) SUM_DEN=SUM_DEN+CI_ZM1_PI_INV ENDDO IF(EPS.NE.ZERO) THEN RES=EXPM1(Z_M_0P5*LOG1P(EPS/Z_PG_M0P5) & +EPS*LOG(EPS_PZ_PG_M0P5)-EPS+LOG1P(-EPS*SUM_NUM/SUM_DEN))& /EPS RETURN ELSE RES=Z_M_0P5/Z_PG_M0P5 & +LOG(EPS_PZ_PG_M0P5)-ONE-SUM_NUM/SUM_DEN RETURN ENDIF ELSE IF(EPS.NE.ZERO) THEN PI_EPS=M_PI*EPS TT=SIN(PI_EPS)/TANZ(M_PI*(Z-N)) T1_EPS_Z=(COS(PI_EPS)+TT)*& GAMMA_RATIO_DIFF_SMALL_EPS(ONE-Z,-EPS) SIN_PI_2_EPS=SIN(M_PI_2*EPS) T2_EPS_Z=(TWO*SIN_PI_2_EPS*SIN_PI_2_EPS-TT)/EPS T_EPS_Z=T1_EPS_Z+T2_EPS_Z RES=(T_EPS_Z/(ONE-EPS*T_EPS_Z)) RETURN ELSE RES=GAMMA_RATIO_DIFF_SMALL_EPS(ONE-Z,-EPS) & -M_PI/TANZ(M_PI*(Z-N)) RETURN ENDIF ENDIF END FUNCTION GAMMA_RATIO_DIFF_SMALL_EPS FUNCTION GAMMA_INV_DIFF_EPS(Z,EPS) ! Calculation of G(z,eps) = [Gamma_inv(z) - Gamma_inv(z+eps)]/eps ! with e and z complex ! ! The G(z,eps) function is used in formulas occuring in 1-z ! and 1/z transformations (see Comp. Phys. Comm. paper). ! Several case have to be considered for its evaluation. ! eps is considered equal to zero ! if z+eps and z are equal numerically. ! ! |eps|oo > 0.1 ! ------------- ! A direct evaluation with the values Gamma_inv(z) ! and Gamma_inv(z+eps) is stable and returned. ! ! |eps|oo <= 0.1 with z+eps and z numerically different ! ----------------------------------------------------- ! If z is a negative integer, z+eps is not, ! so that G(z,eps) = -Gamma_inv(z+eps)/eps, ! for which a direct evaluation is precise and returned. ! If z+eps is a negative integer, z is not, ! so that G(z,eps) = Gamma_inv(z)/eps, ! for which a direct evaluation is precise and returned. ! If both of them are not negative integers, ! one looks for the one of z and z+eps ! which is the closest to a negative integer. ! If it is z, one returns H(z,eps).Gamma_inv(z+eps). ! If it is z+eps, one returns H(z+eps,-eps).Gamma_inv(z). ! Both values are equal, so that one chooses the one ! which makes the Gamma ratio Gamma(z+eps)/Gamma(z) ! in H(z,eps) the smallest in modulus. ! ! z+eps and z numerically equal ! ----------------------------- ! If z is negative integer, G(z,0) = (-1)^(n+1) n!, ! where z = -n, n integer, which is returned. ! If z is not negative integer, one returns H(z,eps).Gamma_inv(z+eps) ! ! Variables ! --------- ! z,eps: input variables of the function G(z,eps) ! eps_pz,x,eps_px: z+eps,real parts of z and z+eps. ! n,m: closest integer ot the real part of z, same for z+eps. ! fact,k: (-1)^(n+1) n!, returned when z = -n, n integer ! and z and z+eps identical numerically (eps ~ 0). ! It is calculated with integer index k. ! is_z_negative_integer,is_eps_pz_negative_integer: ! true if z is a negative integer, false if not, same for z+eps. ! z_neg_int_distance, eps_pz_neg_int_distance: ! |z + |n||oo, |z + eps + |m||oo. ! If |z + |n||oo < |z + eps + |m||oo, ! z is closer to the set of negative integers than z+eps. ! Gamma_inv(z+eps) is then of moderate modulus ! if Gamma_inv(z) is very small. ! If z ~ n, H(z,eps) ~ -1/eps, ! that so returning ! G(z,eps) = H(z,eps).Gamma_inv(z+eps) here is preferred. ! Same for |z + |n||oo > |z + eps + |m||oo with z <-> z+eps. !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: Z,EPS INTEGER(i4) :: M,N,K REAL(wp) :: X,EPS_PX,FACT REAL(wp) :: Z_NEG_INT_DISTANCE REAL(wp) :: EPS_PZ_NEG_INT_DISTANCE COMPLEX(wp) :: GAMMA_INV_DIFF_EPS,EPS_PZ LOGICAL :: IS_Z_NEG_INT,IS_EPS_PZ_NEG_INT EPS_PZ=Z+EPS; X=REAL(Z,wp); EPS_PX=REAL(EPS_PZ,wp) N=NINT(X); M=NINT(EPS_PX) IS_Z_NEG_INT=(Z.EQ.N).AND.(N.LE.0) IS_EPS_PZ_NEG_INT=(EPS_PZ.EQ.M).AND.(M.LE.0) IF(INF_NORM(EPS).GT.0.10_wp) THEN GAMMA_INV_DIFF_EPS = (GAMMA_INV (Z) - GAMMA_INV (EPS_PZ))/EPS RETURN ELSE IF(EPS_PZ.NE.Z) THEN IF(IS_Z_NEG_INT) THEN GAMMA_INV_DIFF_EPS = (-GAMMA_INV (EPS_PZ)/EPS) RETURN ELSE IF(IS_EPS_PZ_NEG_INT) THEN GAMMA_INV_DIFF_EPS = (GAMMA_INV (Z)/EPS) RETURN ELSE Z_NEG_INT_DISTANCE = INF_NORM (Z + ABS (N)) EPS_PZ_NEG_INT_DISTANCE = INF_NORM (EPS_PZ + ABS (M)) IF(Z_NEG_INT_DISTANCE.LT.EPS_PZ_NEG_INT_DISTANCE) THEN GAMMA_INV_DIFF_EPS= & GAMMA_RATIO_DIFF_SMALL_EPS (Z,EPS)*GAMMA_INV (EPS_PZ) RETURN ELSE GAMMA_INV_DIFF_EPS= & GAMMA_RATIO_DIFF_SMALL_EPS (EPS_PZ,-EPS)*GAMMA_INV (Z) RETURN ENDIF ENDIF ELSE IF(IS_Z_NEG_INT.AND.IS_EPS_PZ_NEG_INT) THEN FACT = -ONE;K=-1 DO WHILE (K.GE.N) FACT=FACT*K K=K-1 ENDDO GAMMA_INV_DIFF_EPS = FACT RETURN ELSE GAMMA_INV_DIFF_EPS = & GAMMA_RATIO_DIFF_SMALL_EPS (Z,EPS)*GAMMA_INV (EPS_PZ) RETURN ENDIF END FUNCTION GAMMA_INV_DIFF_EPS FUNCTION A_SUM_INIT(M,EPS,GAMMA_INV_ONE_MEPS) ! Calculation of Gamma_inv(1-m-eps)/eps of the A(z) polynomial in 1-z ! and 1/z transformations ! ! This value occurs in A(z) in 1-z and 1/z transformations ! (see Comp. Phys. Comm. paper) for m > 0. ! Both cases of 1-m-eps numerically negative integer ! or not have to be considered ! ! 1-eps-m and 1-m numerically different ! ------------------------------------- ! One returns Gamma_inv(1-m-eps)/eps directly ! as its value is accurate. ! To calculate Gamma_inv(1-m-eps), ! one uses the value Gamma_inv(1-eps), ! needed in considered transformations, ! and one uses the equality ! Gamma_inv(1-m-eps) = Gamma_inv(1-eps) \prod_{i=1}^{m} (1-eps-i) ! for m > 0. ! It is trivially demonstrated ! from the equality Gamma(x+1) = x.Gamma(x). ! One Gamma function evaluation is removed this way ! from the calculation. ! ! 1-eps-m and 1-m numerically equal ! --------------------------------- ! This implies that 1-m-eps is negative integer numerically. ! Here, eps~0, so that one returns the limit of Gamma_inv(1-m-eps)/eps ! for eps -> 0, which is (-1)^m (m-1)! ! ! Variables ! --------- ! m,eps: variable inputs of the function ! (m,eps) -> Gamma_inv(1-m-eps)/eps ! Gamma_inv_one_meps: Gamma_inv(1-eps), ! previously calculated and here recycled ! to quickly calculate Gamma_inv(1-m-eps). ! one_meps: 1-eps !---------------------------------------------------------------------- INTEGER(i4),INTENT(IN) :: M COMPLEX(wp),INTENT(IN) :: EPS,GAMMA_INV_ONE_MEPS INTEGER(i4) :: N,I REAL(wp) :: FACT COMPLEX(wp) :: A_SUM_INIT,ONE_MEPS COMPLEX(wp) :: GAMMA_INV_ONE_MEPS_MM ! ONE_MEPS = ONE - EPS IF(ONE_MEPS-M.NE.1-M) THEN GAMMA_INV_ONE_MEPS_MM = GAMMA_INV_ONE_MEPS DO I=1,M GAMMA_INV_ONE_MEPS_MM = GAMMA_INV_ONE_MEPS_MM*(ONE_MEPS-I) ENDDO A_SUM_INIT=GAMMA_INV_ONE_MEPS_MM/EPS RETURN ELSE FACT=ONE DO N=2,M-1 FACT=FACT*N ENDDO IF(MOD(M,2).EQ.0) THEN A_SUM_INIT=FACT ELSE A_SUM_INIT=-FACT ENDIF RETURN ENDIF END FUNCTION A_SUM_INIT FUNCTION LOG_A_SUM_INIT(M,EPS) ! Calculation of the log of Gamma_inv(1-m-eps)/eps ! ! See previous function. ! It is used in case Gamma_inv(1-m-eps)/eps might overflow. ! ! Variables ! --------- ! m,eps: variable inputs of the function ! (m,eps) -> log[Gamma_inv(1-m-eps)/eps] ! one_meps_mm: 1-eps-m ! i_Pi: i.Pi ! log_fact: logarithm of (-1)^m (m-1)!, ! here defined as log((m-1)!) + i.Pi if m is odd. !---------------------------------------------------------------------- INTEGER(i4),INTENT(IN) :: M COMPLEX(wp),INTENT(IN) :: EPS INTEGER(i4) :: N REAL(wp) :: LOG_FACT COMPLEX(wp) :: ONE_MEPS_MM,LOG_A_SUM_INIT ONE_MEPS_MM=ONE-EPS-M IF(ONE_MEPS_MM.NE.1-M) THEN LOG_A_SUM_INIT=(-LOG_GAMMA(ONE_MEPS_MM) - LOG(EPS)) RETURN ELSE LOG_FACT=ZERO DO N=2,M-1 LOG_FACT=LOG_FACT + LOG(DBLE(N)) ENDDO IF(MOD(M,2).EQ.0) THEN LOG_A_SUM_INIT=LOG_FACT ELSE LOG_A_SUM_INIT=CMPLX(LOG_FACT,M_PI,wp) ENDIF RETURN ENDIF END FUNCTION LOG_A_SUM_INIT FUNCTION B_SUM_INIT_PS_ONE( & A,B,GAMMA_C,GAMMA_INV_ONE_MEPS,GAMMA_INV_EPS_PA_PM,GAMMA_INV_EPS_PB_PM,MZP1,M,EPS & ) ! Calculation of the first term of the B(z) power series ! in the 1-z transformation, divided by (1-z)^m ! ! In the 1-z transformation, ! the power series B(z) = \sum_{n=0}^{+oo} \beta_n (1-z)^n occurs ! (see Comp. Phys. Comm. paper). ! The first term \beta_0, divided by (1-z)^m, is calculated here. ! m is the closest integer to Re(c-a-b) >= 0 and eps = c-a-b-m. ! ! One has to consider |eps|oo > 0.1 and |eps|oo <= 0.1, ! where 1-m-eps and 1-m can be different or equal numerically, ! leading to some changes in this last case. ! ! |eps|oo > 0.1 ! ------------- ! One has \beta_0/(1-z)^m = [(a)_m (b)_m Gamma_inv(1-eps) ! Gamma_inv(a+m+eps) Gamma_inv(b+m+eps) Gamma_inv(m+1) ! - (1-z)^eps Gamma_inv(a) Gamma_inv(b) Gamma_inv(1+m+eps)] ! [Gamma(c)/eps], stable in this regime for a direct evaluation. ! ! The values of Gamma(c), Gamma_inv(a+m+eps) ! and Gamma_inv(b+m+eps) were already calculated and recycled here. ! Gamma_inv(m+1) is calculated as 1/(m!). ! ! Gamma_inv(1+m+eps) is calculated from Gamma_inv(1-eps), ! using the equalities: ! Gamma_inv(1-m-eps) = Gamma_inv(1-eps) \prod_{i=1}^{m} (1-eps-i), ! where the product is 1 by definition if m = 0, ! Gamma_inv(1+m+eps) = (-1)^m sin (pi.eps) ! /[pi.(eps+m).Gamma_inv(1-m-eps)] ! from Euler reflection formula, Gamma(x+1) = x.Gamma(x) equality, ! and m+eps no zero. ! This scheme is much faster than ! to recalculate Gamma_inv(1+m+eps) directly. ! ! |eps|oo <= 0.1 ! -------------- ! The \beta_0/(1-z)^m expression is rewritten ! so that it contains no instabilities: ! \beta_0/(1-z)^m = Gamma_inv(a+m+eps) Gamma_inv(b+m+eps) ! [(G(1,-eps) Gamma_inv(m+1) + G(m+1,eps)) ! - Gamma_inv(1+m+eps) (G(a+m,eps) Gamma_inv(b+m+eps) ! + G(b+m,eps) Gamma_inv(a+m)) ! - E(log(1-z),eps) Gamma_inv(a+m) Gamma_inv(b+m) Gamma_inv(1+m+eps)] ! (a)_m (b)_m Gamma(c) ! ! E(log(1-z),eps) is [(1-z)^eps - 1]/eps ! if 1-m-eps and 1-m are different numerically, ! and log(1-z) otherwise (eps ~ 0). ! If 1-m-eps and 1-m are equal numerically, ! Gamma_inv(1+m+eps) is numerically equal to Gamma_inv(1+m), ! already calculated as 1/(m!). ! See |eps|oo > 0.1 case for data recycling of other values ! or for 1-m-eps and 1-m different numerically. ! !---------------------------------------------------------------------- ! Variables ! --------- ! a,b,c,one_minus_z: a,b,c and 1-z parameters and arguments ! of the 2F1(a,b,c,z) function. ! m,eps: closest integer to c-a-b, with Re(c-a-b) >= 0 ! and eps = c-a-b-m ! Gamma_c,Gamma_inv_one_meps,Gamma_inv_eps_pa_pm, Gamma_inv_eps_pb_pm: ! recycled values of Gamma(c), Gamma_inv(1-eps), ! Gamma_inv(a+m+eps) and Gamma_inv(b+m+eps). ! inf_norm_eps,phase,a_pm,b_pm,one_meps,Pi_eps,Pi_eps_pm: ! |eps|oo,(-1)^m,a+m,b+m,1-eps,pi.eps,pi.(eps+m) ! Gamma_inv_one_meps_mm,Gamma_inv_eps_pm_p1: ! Gamma_inv(1-m-eps) and Gamma_inv(1+m+eps) ! calculated with the recycling scheme. ! prod1: (a)_m (b)_m Gamma_inv(1-eps) Gamma_inv(a+m+eps) ! x Gamma_inv(b+m+eps) Gamma_inv(m+1) in |eps|oo > 0.1 case. ! prod2: (1-z)^eps Gamma_inv(a) Gamma_inv(b) Gamma_inv(1+m+eps) ! in |eps|oo > 0.1 case. ! Gamma_inv_mp1,prod_ab: Gamma_inv(m+1) calculated as 1/(m!) ! and (a)_m (b)_m in |eps|oo <= 0.1 case. ! is_eps_non_zero: true if 1-m-eps and 1-m are different numerically, ! false if not. ! Gamma_inv_a_pm,Gamma_inv_b_pm,z_term: Gamma_inv(a+m),Gamma_inv(b+m), ! E(eps,log(1-z)) ! prod1: Gamma_inv(a+m+eps) Gamma_inv(b+m+eps) ! x [(G(1,-eps) Gamma_inv(m+1) + G(m+1,eps)) in |eps|oo <= 0.1 case. ! prod2: Gamma_inv(1+m+eps) (G(a+m,eps) Gamma_inv(b+m+eps) ! + G(b+m,eps) Gamma_inv(a+m)) ! prod3: E(eps,log(1-z)) Gamma_inv(a+m) Gamma_inv(b+m) ! Gamma_inv(1+m+eps) ! res: returned \beta_0/(1-z)^m value in all cases. !---------------------------------------------------------------------- INTEGER(i4),INTENT(IN) :: M COMPLEX(wp),INTENT(IN) :: A,B,GAMMA_C,GAMMA_INV_ONE_MEPS, & GAMMA_INV_EPS_PA_PM,GAMMA_INV_EPS_PB_PM,MZP1,EPS INTEGER(i4) :: M_M1,N,I,PHASE REAL(wp) :: INF_NORM_EPS,GAMMA_INV_MP1 COMPLEX(wp) :: A_PM,B_SUM_INIT_PS_ONE,PI_EPS,GAMMA_INV_ONE_MEPS_MM COMPLEX(wp) :: B_PM,TMP1,TMP2 COMPLEX(wp) :: Z_TERM,PROD1,PROD2,PROD3,ONE_MEPS,PI_EPS_PM COMPLEX(wp) :: GAMMA_INV_A_PM,PROD_AB,GAMMA_INV_B_PM COMPLEX(wp) :: GAMMA_INV_EPS_PM_P1 INF_NORM_EPS=INF_NORM(EPS); M_M1=M-1; A_PM=A+M; B_PM=B+M ONE_MEPS=ONE-EPS; PI_EPS=M_PI*EPS; PI_EPS_PM = M_PI*(EPS+M) IF(MOD(M,2).EQ.0) THEN PHASE = 1 ELSE PHASE = -1 ENDIF GAMMA_INV_ONE_MEPS_MM = GAMMA_INV_ONE_MEPS DO I=1,M GAMMA_INV_ONE_MEPS_MM = GAMMA_INV_ONE_MEPS_MM*(ONE_MEPS - I) ENDDO IF(INF_NORM_EPS.GT.0.10_wp) THEN GAMMA_INV_EPS_PM_P1 = PHASE*SIN(PI_EPS) & /(PI_EPS_PM*GAMMA_INV_ONE_MEPS_MM) PROD1=GAMMA_INV_ONE_MEPS*GAMMA_INV_EPS_PA_PM*GAMMA_INV_EPS_PB_PM DO N=0,M_M1 PROD1=PROD1*(A+N)*(B+N)/(N+ONE) ENDDO PROD2=GAMMA_INV(A)*GAMMA_INV(B)*GAMMA_INV_EPS_PM_P1*(MZP1**EPS) B_SUM_INIT_PS_ONE=GAMMA_C*(PROD1-PROD2)/EPS RETURN ELSE GAMMA_INV_MP1=ONE;PROD_AB=ONE DO N=0,M_M1 GAMMA_INV_MP1 = GAMMA_INV_MP1/(N+ONE) PROD_AB = PROD_AB*(A+N)*(B+N) ENDDO IF(ONE_MEPS-M.NE.1-M) THEN Z_TERM=EXPM1(EPS*LOG(MZP1))/EPS GAMMA_INV_EPS_PM_P1 = PHASE*SIN(PI_EPS) & /(PI_EPS_PM*GAMMA_INV_ONE_MEPS_MM) ELSE Z_TERM=LOG(MZP1) GAMMA_INV_EPS_PM_P1 = GAMMA_INV_MP1 ENDIF GAMMA_INV_A_PM=GAMMA_INV(A_PM);GAMMA_INV_B_PM=GAMMA_INV(B_PM) TMP1=ONE; TMP2=M+1; PROD1 = GAMMA_INV_EPS_PA_PM*GAMMA_INV_EPS_PB_PM & *(GAMMA_INV_MP1*GAMMA_INV_DIFF_EPS(TMP1,-EPS) & +GAMMA_INV_DIFF_EPS(TMP2,EPS)) PROD2 = GAMMA_INV_EPS_PM_P1 & *(GAMMA_INV_EPS_PB_PM*GAMMA_INV_DIFF_EPS(A_PM,EPS) & +GAMMA_INV_A_PM*GAMMA_INV_DIFF_EPS (B_PM,EPS)) PROD3 = GAMMA_INV_A_PM*GAMMA_INV_B_PM*GAMMA_INV_EPS_PM_P1*Z_TERM B_SUM_INIT_PS_ONE=GAMMA_C*PROD_AB*(PROD1-PROD2-PROD3) RETURN ENDIF END FUNCTION B_SUM_INIT_PS_ONE FUNCTION B_SUM_INIT_PS_INFINITY( & A,C,GAMMA_C,GAMMA_INV_CMA,GAMMA_INV_ONE_MEPS,GAMMA_INV_EPS_PA_PM,Z,M,EPS & ) ! Calculation of the first term of the B(z) power series ! in the 1/z transformation, divided by z^{-m} ! ! In the 1/z transformation, the power series ! B(z) = \sum_{n=0}^{+oo} \beta_n z^{-n} occurs ! (see Comp. Phys. Comm. paper). ! The first term \beta_0, divided by z^{-m}, is calculated here. ! m is the closest integer to Re(b-a) >= 0 and eps = b-a-m. ! ! One has to consider |eps|oo > 0.1 and |eps|oo <= 0.1, ! where 1-m-eps and 1-m can be different or equal numerically, ! leading to some changes in this last case. ! ! |eps|oo > 0.1 ! ------------- ! One has \beta_0/z^{-m} = [(a)_m (1-c+a)_m Gamma_inv(1-eps) ! Gamma_inv(a+m+eps) Gamma_inv(c-a) Gamma_inv(m+1) ! - (-z)^{-eps} (1-c+a+eps)_m Gamma_inv(a) Gamma_inv(c-a-eps) ! Gamma_inv(1+m+eps)].[Gamma(c)/eps], ! stable in this regime for a direct evaluation. ! ! The values of Gamma(c), Gamma_inv(c-a) and Gamma_inv(a+m+eps) ! were already calculated and recycled here. ! Gamma_inv(m+1) is calculated as 1/(m!). ! Gamma_inv(1+m+eps) is calculated from Gamma_inv(1-eps) ! as in the 1-z transformation routine. ! ! |eps|oo <= 0.1 ! -------------- ! The \beta_0/z^{-m} expression is rewritten ! so that it contains no instabilities: ! \beta_0/z^{-m} = [((1-c+a+eps)_m G(1,-eps) - P(m,eps,1-c+a) ! Gamma_inv(1-eps)) Gamma_inv(c-a) Gamma_inv(a+m+eps) Gamma_inv(m+1) ! + (1-c+a+eps)_m [G(m+1,eps) Gamma_inv(c-a) Gamma_inv(a+m+eps) ! - G(a+m,eps) Gamma_inv(c-a) Gamma_inv(m+1+eps)] ! - (G(c-a,-eps) - E(log(-z),-eps)) Gamma_inv(m+1+eps) ! Gamma_inv(a+m)]] (a)_m Gamma(c) ! ! Definitions and method are the same ! as in the 1-z transformation routine, except for P(m,eps,1-c+a). ! P(m,eps,s) = [(s+eps)_m - (s)_m]/eps ! for eps non zero and has a limit for eps -> 0. ! Let n0 be the closest integer to -Re(s) for s complex. ! A stable formula available for eps -> 0 for P(m,eps,s) is: ! P(m,eps,s) = (s)_m E(\sum_{n=0}^{m-1} L(1/(s+n),eps),eps) ! if n0 is not in [0:m-1], ! P(m,eps,s) = \prod_{n=0, n not equal to n0}^{m-1} (s+eps+n) ! + (s)_m E(\sum_{n=0, n not equal to n0}^{m-1} L(1/(s+n),eps),eps) ! if n0 is in [0:m-1]. ! L(s,eps) is log1p(s eps)/eps if eps is not zero, ! and L(s,0) = s. ! This expression is used in the code. ! ! Variables ! --------- ! a,b,c,z: a,b,c and z parameters ! and arguments of the 2F1(a,b,c,z) function. ! m,eps: closest integer to b-a, with Re(b-a) >= 0 and eps = b-a-m. ! Gamma_c,Gamma_inv_cma,Gamma_inv_one_meps,Gamma_inv_eps_pa_pm: ! recycled values of Gamma(c), Gamma_inv(c-a), Gamma_inv(1-eps) ! and Gamma_inv(a+m+eps). ! inf_norm_eps,phase,cma,a_mc_p1,a_mc_p1_pm,cma_eps,eps_pa_mc_p1,a_pm: ! |eps|oo,(-1)^m,c-a,1-c+a+m,c-a-eps,1-c+a+eps,a+m ! Gamma_inv_cma_meps,one_meps,Pi_eps,Pi_eps_pm: ! Gamma_inv(c-a-eps),1-eps,pi.eps,pi.(eps+m) ! Gamma_inv_one_meps_mm,Gamma_inv_eps_pm_p1: Gamma_inv(1-m-eps) ! and Gamma_inv(1+m+eps) calculated with the recycling scheme. ! prod1: (a)_m (1-c+a)_m Gamma_inv(1-eps) Gamma_inv(a+m+eps) ! x Gamma_inv(c-a) Gamma_inv(m+1) in |eps|oo > 0.1 case. ! prod2: (-z)^{-eps} (1-c+a+eps)_m Gamma_inv(a) ! x Gamma_inv(c-a-eps) Gamma_inv(1+m+eps) in |eps|oo > 0.1 case. ! n0: closest integer to -Re(1-c+a) ! is_n0_here: true is n0 belongs to [0:m-1], false if not. ! is_eps_non_zero: true if 1-m-eps and 1-m are different numerically, ! false if not. ! Gamma_inv_mp1,prod_a,prod_a_mc_p1: ! Gamma_inv(m+1) calculated as 1/(m!), ! (a)_m and (1-c+a)_m in |eps|oo <= 0.1 case. ! prod_eps_pa_mc_p1_n0: ! \prod_{n=0, n not equal to n0}^{m-1} (1-c+a+eps+n) ! if n0 belongs to [0:m-1], 0.0 if not, in |eps|oo <= 0.1 case. ! prod_eps_pa_mc_p1: (1-c+a+eps)_m in |eps|oo <= 0.1 case. ! sum: \sum_{n=0, n not equal to n0}^{m-1} L(1/(s+n),eps) if 1-m-eps ! and 1-m are different numerically, ! \sum_{n=0, n not equal to n0}^{m-1} 1/(s+n) if not. ! a_pn,a_mc_p1_pn,eps_pa_mc_p1_pn: a+n,1-c+a+n,1-c+a+eps+n values ! used in (a)_m, (1-c+a)_m and (1-c+a+eps)_m evaluations. ! sum_term,prod_diff_eps,z_term: ! E(\sum_{n=0, n not equal to n0}^{m-1} L(1/(s+n),eps),eps), ! P(m,eps,1-c+a), -E(-eps,log(-z)) ! Gamma_inv_a_pm,Gamma_prod1: Gamma_inv(a+m), ! Gamma_inv(c-a).Gamma_inv(a+m+eps) ! prod1: ((1-c+a+eps)_m G(1,-eps) ! - P(m,eps,1-c+a) Gamma_inv(1-eps)) Gamma_inv(c-a) ! x Gamma_inv(a+m+eps) Gamma_inv(m+1) ! prod_2a: Gamma_inv(c-a).Gamma_inv(a+m+eps).G(m+1,eps) ! prod_2b: G(a+m,eps) Gamma_inv(c-a) Gamma_inv(m+1+eps) ! prod_2c: (G(c-a,-eps) ! - E(log(-z),-eps)) Gamma_inv(m+1+eps) Gamma_inv(a+m) ! prod2: (1-c+a+eps)_m [G(m+1,eps) Gamma_inv(c-a) Gamma_inv(a+m+eps) ! - G(a+m,eps) Gamma_inv(c-a) Gamma_inv(m+1+eps)] ! - (G(c-a,-eps) - E(log(-z),-eps)) ! x Gamma_inv(m+1+eps) Gamma_inv(a+m)]] ! res: returned \beta_0/z^{-m} value in all cases. !---------------------------------------------------------------------- INTEGER(i4),INTENT(IN) :: M COMPLEX(wp),INTENT(IN) :: A,C,GAMMA_C,GAMMA_INV_CMA,Z,EPS COMPLEX(wp),INTENT(IN) :: GAMMA_INV_ONE_MEPS,GAMMA_INV_EPS_PA_PM INTEGER(i4) :: M_M1,I,N,N0,PHASE LOGICAL :: IS_N0_HERE,IS_EPS_NON_ZERO REAL(wp) :: INF_NORM_EPS,NP1,GAMMA_INV_MP1 COMPLEX(wp) :: B_SUM_INIT_PS_INFINITY,TMP1 COMPLEX(wp) :: CMA,A_MC_P1,A_MC_P1_PM,CMA_MEPS,EPS_PA_MC_P1,A_PM COMPLEX(wp) :: GAMMA_INV_EPS_PM_P1,GAMMA_INV_CMA_MEPS,PI_EPS COMPLEX(wp) :: PROD1,PROD2,A_PN,A_MC_P1_PN,ONE_MEPS COMPLEX(wp) :: PROD_A,PROD_A_MC_P1,PROD_EPS_PA_MC_P1_N0,PI_EPS_PM COMPLEX(wp) :: PROD_EPS_PA_MC_P1,SUM_N0,Z_TERM,SUM_TERM COMPLEX(wp) :: PROD_DIFF_EPS,GAMMA_INV_A_PM,GAMMA_PROD1 COMPLEX(wp) :: PROD_2A,PROD_2B,PROD_2C COMPLEX(wp) :: EPS_PA_MC_P1_PN,GAMMA_INV_ONE_MEPS_MM ! INF_NORM_EPS=INF_NORM(EPS); CMA=C-A; A_MC_P1=A-C+ONE A_MC_P1_PM=A_MC_P1+M; CMA_MEPS=CMA-EPS; EPS_PA_MC_P1=EPS+A_MC_P1 A_PM=A+M; M_M1=M-1; ONE_MEPS=ONE-EPS; PI_EPS=M_PI*EPS PI_EPS_PM=M_PI*(EPS+M); GAMMA_INV_CMA_MEPS=GAMMA_INV(CMA_MEPS) IF(MOD(M,2).EQ.0) THEN PHASE = 1 ELSE PHASE = -1 ENDIF GAMMA_INV_ONE_MEPS_MM = GAMMA_INV_ONE_MEPS DO I=1,M GAMMA_INV_ONE_MEPS_MM = GAMMA_INV_ONE_MEPS_MM*(ONE_MEPS - I) ENDDO IF(INF_NORM_EPS.GT.0.1_wp) THEN GAMMA_INV_EPS_PM_P1 = PHASE*SIN(PI_EPS) & /(PI_EPS_PM*GAMMA_INV_ONE_MEPS_MM) PROD1 = GAMMA_INV_CMA*GAMMA_INV_EPS_PA_PM*GAMMA_INV_ONE_MEPS PROD2 = GAMMA_INV(A)*GAMMA_INV_CMA_MEPS*GAMMA_INV_EPS_PM_P1 & *((-Z)**(-EPS)) DO N=0,M_M1 A_PN=A+N; A_MC_P1_PN=A_MC_P1+N EPS_PA_MC_P1_PN=EPS+A_MC_P1_PN;NP1=N+ONE PROD1 = PROD1*A_PN*A_MC_P1_PN/NP1 PROD2 = PROD2*EPS_PA_MC_P1_PN ENDDO B_SUM_INIT_PS_INFINITY = GAMMA_C*(PROD1-PROD2)/EPS RETURN ELSE N0=-NINT(REAL(A_MC_P1,wp)) IS_EPS_NON_ZERO=ONE_MEPS-M.NE.1-M IS_N0_HERE=(N0.GE.0).AND.(N0.LT.M) GAMMA_INV_MP1=ONE; PROD_A=ONE; PROD_A_MC_P1=ONE PROD_EPS_PA_MC_P1=ONE; SUM_N0=ZERO IF(IS_N0_HERE) THEN PROD_EPS_PA_MC_P1_N0 = ONE ELSE PROD_EPS_PA_MC_P1_N0 = ZERO ENDIF DO N=0,M_M1 A_PN=A+N; A_MC_P1_PN=A_MC_P1+N EPS_PA_MC_P1_PN=EPS+A_MC_P1_PN; NP1=N+ONE PROD_A = PROD_A*A_PN PROD_A_MC_P1 = PROD_A_MC_P1*A_MC_P1_PN PROD_EPS_PA_MC_P1 = PROD_EPS_PA_MC_P1*EPS_PA_MC_P1_PN GAMMA_INV_MP1 = GAMMA_INV_MP1/NP1 IF(N.NE.N0) THEN IF(IS_N0_HERE) THEN PROD_EPS_PA_MC_P1_N0=PROD_EPS_PA_MC_P1_N0 & *EPS_PA_MC_P1_PN ENDIF IF(IS_EPS_NON_ZERO) THEN SUM_N0 = SUM_N0 + LOG1P(EPS/A_MC_P1_PN) ELSE SUM_N0 = SUM_N0 + ONE/A_MC_P1_PN ENDIF ENDIF ENDDO IF(IS_EPS_NON_ZERO) THEN GAMMA_INV_EPS_PM_P1 = PHASE*SIN(PI_EPS) & /(PI_EPS_PM*GAMMA_INV_ONE_MEPS_MM) SUM_TERM = EXPM1(SUM_N0)/EPS Z_TERM = EXPM1(-EPS*LOG(-Z))/EPS ELSE GAMMA_INV_EPS_PM_P1 = GAMMA_INV_MP1 SUM_TERM = SUM_N0 Z_TERM = -LOG(-Z) ENDIF PROD_DIFF_EPS = PROD_EPS_PA_MC_P1_N0 + PROD_A_MC_P1*SUM_TERM GAMMA_INV_A_PM = GAMMA_INV(A_PM) GAMMA_PROD1=GAMMA_INV_CMA*GAMMA_INV_EPS_PA_PM TMP1=ONE PROD1 = GAMMA_PROD1*GAMMA_INV_MP1*(GAMMA_INV_DIFF_EPS(TMP1,-EPS) & *PROD_EPS_PA_MC_P1 - GAMMA_INV_ONE_MEPS*PROD_DIFF_EPS) TMP1=M+1 PROD_2A = GAMMA_PROD1*GAMMA_INV_DIFF_EPS(TMP1,EPS) PROD_2B = GAMMA_INV_CMA*GAMMA_INV_EPS_PM_P1 & *GAMMA_INV_DIFF_EPS(A_PM,EPS) PROD_2C = GAMMA_INV_EPS_PM_P1*GAMMA_INV_A_PM & *(GAMMA_INV_DIFF_EPS(CMA,-EPS) + GAMMA_INV_CMA_MEPS*Z_TERM) PROD2 = PROD_EPS_PA_MC_P1*(PROD_2A - PROD_2B - PROD_2C) B_SUM_INIT_PS_INFINITY = GAMMA_C*PROD_A*(PROD1+PROD2) RETURN ENDIF END FUNCTION B_SUM_INIT_PS_INFINITY SUBROUTINE CV_POLY_DER_TAB_CALC(A,B,C,Z,CV_POLY_DER_TAB) ! Calculation of the derivative of the polynomial P(X) ! ! testing power series convergence ! ! P(X) = |z(a+X)(b+X)|^2 - |(c+X)(X+1)|^2 ! = \sum_{i=0}^{4} c[i] X^{i}, for |z| < 1. ! It is positive when the power series term modulus increases ! and negative when it decreases, ! so that its derivative provides information on its convergence ! (see Comp. Phys. Comm. paper). ! Its derivative components cv_poly_der_tab[i] = (i+1) c[i+1] ! for i in [0:3] ! so that P'(X) = \sum_{i=0}^{3} cv_poly_der_tab[i] X^{i} ! are calculated. ! ! Variables: ! ---------- ! a,b,c,z: a,b,c and z parameters and arguments ! of the 2F1(a,b,c,z) function. ! cv_poly_der_tab[3]: table of four doubles ! containing the P'(X) components. ! mod_a2,mod_b2,mod_c2,mod_z2,R_a,Re_b,Re_c: |a|^2, |b|^2, |c|^2, ! |z|^2, Re(a), Re(b), Re(c), with which P(X) can be expressed. !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: A,B,C,Z REAL(wp),INTENT(OUT) :: CV_POLY_DER_TAB(0:3) REAL(wp) :: MOD_A2,MOD_B2,MOD_C2,MOD_Z2 REAL(wp) :: RE_A,RE_B,RE_C,IM_A,IM_B,IM_C,RE_Z,IM_Z ! RE_A=REAL(A,wp); IM_A=AIMAG(A); MOD_A2=RE_A*RE_A+IM_A*IM_A RE_B=REAL(B,wp); IM_B=AIMAG(B); MOD_B2=RE_B*RE_B+IM_B*IM_B RE_C=REAL(C,wp); IM_C=AIMAG(C); MOD_C2=RE_C*RE_C+IM_C*IM_C RE_Z=REAL(Z,wp); IM_Z=AIMAG(Z); MOD_Z2=RE_Z*RE_Z+IM_Z*IM_Z CV_POLY_DER_TAB(0)=TWO*((RE_A*MOD_B2+RE_B*MOD_A2)*MOD_Z2-RE_C-MOD_C2) CV_POLY_DER_TAB(1)=TWO*((MOD_A2+MOD_B2+4.0_wp*RE_A*RE_B)*MOD_Z2 & -ONE-4.0_wp*RE_C-MOD_C2) CV_POLY_DER_TAB(2)=6.0_wp*((RE_A+RE_B)*MOD_Z2-RE_C-ONE) CV_POLY_DER_TAB(3)=4.0_wp*(MOD_Z2-ONE) END SUBROUTINE CV_POLY_DER_TAB_CALC FUNCTION CV_POLY_DER_CALC(CV_POLY_DER_TAB,X) ! Calculation of the derivative of the polynomial P(X) ! ! testing power series convergence at one x value ! ! P'(x) is calculated for a real x. ! See P'(X) components calculation routine for definitions. !---------------------------------------------------------------------- REAL(wp),INTENT(IN) :: X REAL(wp),INTENT(IN) :: CV_POLY_DER_TAB(0:3) REAL(wp) :: CV_POLY_DER_CALC ! CV_POLY_DER_CALC=CV_POLY_DER_TAB(0)+X*(CV_POLY_DER_TAB(1) & +X*(CV_POLY_DER_TAB(2)+X*CV_POLY_DER_TAB(3))) RETURN END FUNCTION CV_POLY_DER_CALC FUNCTION MIN_N_CALC(CV_POLY_DER_TAB) ! Calculation of an integer after which false convergence cannot occur ! ! See cv_poly_der_tab_calc routine for definitions. ! If P'(x) < 0 and P''(x) < 0 for x > xc, it will be so for all x > xc ! as P(x) -> -oo for x -> +oo ! and P(x) can have at most one maximum for x > xc. ! It means that the 2F1 power series term modulus will increase ! or decrease to 0 for n > nc, ! with nc the smallest positive integer larger than xc. ! ! If P'(X) = C0 + C1.X + C2.X^2 + C3.X^3, ! the discriminant of P''(X) is Delta = C2^2 - 3 C1 C3. ! ! If Delta > 0, P''(X) has two different real roots ! and its largest root is -(C2 + sqrt(Delta))/(3 C3), ! because C3 = 4(|z|^2 - 1) < 0. ! One can take xc = -(C2 + sqrt(Delta))/(3 C3) ! and one returns its associated nc integer. ! ! If Delta <= 0, P''(X) has at most one real root, ! so that P'(X) has only one root and then P(X) only one maximum. ! In this case, one can choose xc = nc = 0, which is returned. ! ! Variables ! --------- ! cv_poly_der_tab: table of four doubles ! containing the P'(X) coefficients ! C1,C2,three_C3: cv_poly_der_tab[1], cv_poly_der_tab[2] ! and 3.0*cv_poly_der_tab[3], so that P''(X) = C1 + 2.C2.x + three_C3.x^2 ! Delta: discriminant of P''(X), equal to C2^2 - 3 C1 C3. ! largest_root: if Delta > 0, ! P''(X) largest real root equal to -(C2 + sqrt(Delta))/(3 C3). !---------------------------------------------------------------------- REAL(wp),INTENT(IN) :: CV_POLY_DER_TAB(0:3) INTEGER(i4) :: MIN_N_CALC REAL(wp) :: C1,C2,THREE_C3,DELTA,LARGEST_ROOT ! C1=CV_POLY_DER_TAB(1); C2=CV_POLY_DER_TAB(2) THREE_C3=3.0_wp*CV_POLY_DER_TAB(3); DELTA = C2*C2 - THREE_C3*C1 IF(DELTA.LE.ZERO) THEN MIN_N_CALC = 0 RETURN ELSE LARGEST_ROOT = -(C2 + SQRT (DELTA))/THREE_C3 MIN_N_CALC = MAX(CEILING(LARGEST_ROOT),0) RETURN ENDIF END FUNCTION MIN_N_CALC FUNCTION HYP_PS_ZERO(A,B,C,Z) ! Calculation of the 2F1 power series converging for |z| < 1 ! ! One has 2F1(a,b,c,z) ! = \sum_{n = 0}^{+oo} (a)_n (b)_n / ((c)_n n!) z^n, ! so that 2F1(a,b,c,z) = \sum_{n = 0}^{+oo} t[n] z^n, ! with t[0] = 1 and t[n+1] = (a+n)(b+n)/((c+n)(n+1)) t[n] for n >= 0. ! If a or b are negative integers, ! F(z) is a polynomial of degree -a or -b, evaluated directly. ! If not, one uses the test of convergence |t[n] z^n|oo < 1E-15 ! to truncate the series after it was checked ! that false convergence cannot occur. ! ! Variables: ! ---------- ! a,b,c,z: a,b,c and z parameters and arguments ! of the 2F1(a,b,c,z) function. One must have here |z| < 1. ! term,sum: term of the 2F1 power series equal to t[n] z^n, ! truncated sum at given n of the 2F1 power series. ! na,nb: absolute values of the closest integers to Re(a) and Re(b). ! a = -na or b = -nb means one is in the polynomial case. ! cv_poly_der_tab: coefficients of the derivative ! of the polynomial P(X) = |z(a+X)(b+X)|^2 - |(c+X)(X+1)|^2 ! min_n: smallest integer after which false convergence cannot occur. ! It is calculated in min_n_calc. ! possible_false_cv: always true if n < min_n. ! If n >= min_n, it is true if P'(n) > 0. ! If n >= min_n and P'(n) < 0, ! it becomes false and remains as such for the rest of the calculation. ! One can then check if |t[n] z^n|oo < 1E-15 to truncate the series. !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: A,B,C,Z INTEGER(i4) :: N,NA,NB,MIN_N COMPLEX(wp) :: HYP_PS_ZERO,TERM LOGICAL :: POSSIBLE_FALSE_CV REAL(wp) :: CV_POLY_DER_TAB(0:3) NA = ABS(NINT(REAL(A,wp))) NB = ABS(NINT(REAL(B,wp))) TERM=ONE; HYP_PS_ZERO=ONE IF(A.EQ.(-NA)) THEN DO N=0,NA-1 TERM = TERM*Z*(A+N)*(B+N)/((N+ONE)*(C+N)) HYP_PS_ZERO = HYP_PS_ZERO + TERM ENDDO RETURN ELSE IF(B.EQ.(-NB)) THEN DO N=0,NB-1 TERM = TERM*Z*(A+N)*(B+N)/((N+ONE)*(C+N)) HYP_PS_ZERO = HYP_PS_ZERO + TERM ENDDO RETURN ELSE CALL CV_POLY_DER_TAB_CALC(A,B,C,Z,CV_POLY_DER_TAB) POSSIBLE_FALSE_CV=.TRUE. MIN_N=MIN_N_CALC(CV_POLY_DER_TAB);N=0 DO WHILE(POSSIBLE_FALSE_CV.OR.(INF_NORM(TERM).GT.eps_wp)) TERM = TERM*Z*(A+N)*(B+N)/((N+ONE)*(C+N)) HYP_PS_ZERO = HYP_PS_ZERO + TERM IF(POSSIBLE_FALSE_CV.AND.(N.GT.MIN_N)) THEN POSSIBLE_FALSE_CV = & (CV_POLY_DER_CALC(CV_POLY_DER_TAB,DBLE(N)).GT.ZERO) ENDIF N=N+1 ENDDO RETURN ENDIF END FUNCTION HYP_PS_ZERO FUNCTION HYP_PS_ONE(A,B,C,MZP1) ! Calculation of the 2F1 power series ! converging with the 1-z transformation ! ! The formula for F(z) in the 1-z transformation holds: ! F(z) = (-1)^m (pi.eps)/sin (pi.eps) [A(z) + B(z)] ! for eps not equal to zero, F(z) = (-1)^m [A(z) + B(z)] for eps = 0 ! where m = |Re(c-a-b)], eps = c-a-b-m, ! A(z) = \sum_{n=0}^{m-1} alpha[n] (1-z)^n, ! B(z) = \sum_{n=0}^{+oo} beta[n] (1-z)^n, and: ! ! alpha[0] = [Gamma_inv(1-m-eps)/eps] Gamma_inv(a+m+eps) ! x Gamma_inv(b+m+eps) Gamma(c) ! [Gamma_inv(1-m-eps)/eps] is calculated in A_sum_init. ! alpha[0] is calculated with log[Gamma] ! if the previous expression might overflow, ! and its imaginary part removed if a, b and c are real. ! alpha[n+1] = (a+n)(b+n)/[(n+1)(1-m-eps+n)] alpha[n], n in [0:m-2]. ! ! beta[0] is defined in B_sum_init_PS_one function comments. ! gamma[0] = Gamma(c) (a)_m (b)_m (1-z)^m Gamma_inv(a+m+eps) ! x Gamma_inv(b+m+eps) Gamma_inv(m+1) Gamma_inv(1-eps) ! ! beta[n+1] = (a+m+n+eps)(b+m+n+eps)/[(m+n+1+eps)(n+1)] beta[n] ! + [(a+m+n)(b+m+n)/(m+n+1) - (a+m+n) - (b+m+n) - eps ! + (a+m+n+eps)(b+m+n+eps)/(n+1)] ! x gamma[n]/[(n+m+1+eps)(n+1+eps)], n >= 0. ! gamma[n+1] = (a+m+n)(b+m+n)/[(m+n+1)(n+1-eps)] gamma[n], n >= 0. ! ! B(z) converges <=> |1-z| < 1 ! The test of convergence is |beta[n] (1-z)^n|oo < 1E-15 |beta[0]|oo ! for n large enough so that false convergence cannot occur. ! ! Variables ! --------- ! a,b,c,one_minus_z: a,b,c parameters ! and 1-z from z argument of 2F1(a,b,c,z) ! m,phase,m_p1,eps,eps_pm,eps_pm_p1, ! a_pm,b_pm,one_meps,one_meps_pm: ! |Re(c-a-b)], (-1)^m, m+1, c-a-b-m, ! eps+m, eps+m+1, a+m, b+m, 1-eps, 1-eps-m ! eps_pa,eps_pb,eps_pa_pm,eps_pb_pm,Pi_eps,Gamma_c: ! eps+a, eps+b, eps+a+m, eps+b+m, pi.eps, Gamma(c) ! Gamma_inv_eps_pa_pm,Gamma_inv_eps_pb_pm,Gamma_prod: ! Gamma_inv(eps+a+m), Gamma_inv(eps+b+m), ! Gamma(c).Gamma_inv(eps+a+m).Gamma_inv(eps+b+m) ! Gamma_inv_one_meps,A_first_term,A_sum,A_term: ! Gamma_inv(1-eps), alpha[0], A(z), alpha[n] (1-z)^n ! pow_mzp1_m,B_first_term,prod_B,ratio: (1-z)^m, beta[0], ! (a)_m (b)_m (1-z)^m, (a+n)(b+n)/(n+1) for n in [0:m-2]. ! B_extra_term,B_term,B_sum,B_prec: ! gamma[n], beta[n] (1-z)^n, B(z), 1E-15 |beta[0|oo ! cv_poly1_der_tab,cv_poly2_der_tab: P1'(X) and P2'(X) coefficients ! of the potentials derivatives of P1(X) and P2(X) ! defined in cv_poly_der_tab_calc with parameters ! a1 = a, b1 = b, c1 = 1-m-eps, z1 = 1-z ! and a2 = eps+b+m, b2 = eps+a+m,c2 = eps+m+1, z2 = 1-z. ! min_n: smallest integer after which false convergence cannot occur. ! It is calculated in min_n_calc with both P1'(X) and P2'(X), ! so one takes the largest integer coming from both calculations. ! possible_false_cv: always true if n < min_n. ! If n >= min_n, it is true if P1'(n) > 0 or P2'(n) > 0. ! If n >= min_n and P1'(n) < 0 and P2'(n) < 0, ! it becomes false and remains as such for the rest of the calculation. ! One can then check if |beta[n] z^n|oo < 1E-15 to truncate the series. ! n,n_pm_p1,n_p1,a_pm_pn,b_pm_pn,eps_pm_p1_pn,n_p1_meps,eps_pa_pm_pn, ! eps_pb_pm_pn,eps_pm_pn: index of power series, n+m+1, n+1, ! a+m+n, b+m+n, eps+m+n+1, n+1-eps, eps+a+m+n, eps+b+m+n, eps+m+n, ! prod1,prod2,prod3: (eps+a+m+n)(eps+b+m+n), ! (eps+m+1+n)(n+1), (a+m+n)(b+m+n) !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: A,B,C,MZP1 INTEGER(i4) :: N,M,PHASE,M_M2,MIN_N,M_P1 REAL(wp) :: B_PREC,N_P1,N_PM_P1 COMPLEX(wp) :: HYP_PS_ONE,EPS,EPS_PM,EPS_PM_P1,A_PM COMPLEX(wp) :: B_PM,ONE_MEPS_MM,EPS_PA,EPS_PB,PI_EPS,GAMMA_PROD COMPLEX(wp) :: EPS_PA_PM,EPS_PB_PM COMPLEX(wp) :: A_SUM,A_TERM,ONE_MEPS COMPLEX(wp) :: B_EXTRA_TERM,B_TERM,B_SUM,GAMMA_C,RATIO COMPLEX(wp) :: A_PM_PN,B_PM_PN,EPS_PM_P1_PN,N_P1_MEPS COMPLEX(wp) :: PROD1,PROD2,PROD3 COMPLEX(wp) :: EPS_PA_PM_PN,EPS_PB_PM_PN,EPS_PM_PN,PROD_B,POW_MZP1_M COMPLEX(wp) :: GAMMA_INV_EPS_PA_PM,GAMMA_INV_EPS_PB_PM COMPLEX(wp) :: GAMMA_INV_ONE_MEPS LOGICAL :: POSSIBLE_FALSE_CV REAL(wp) :: CV_POLY1_DER_TAB(0:3),CV_POLY2_DER_TAB(0:3) M=NINT(REAL(C-A-B,wp)); M_M2=M-2; M_P1=M+1 IF(MOD(M,2).EQ.0) THEN PHASE=1 ELSE PHASE=-1 ENDIF EPS=C-A-B-M; EPS_PM=EPS+M; EPS_PM_P1=EPS_PM+ONE; A_PM=A+M;B_PM=B+M ONE_MEPS=ONE-EPS; ONE_MEPS_MM=ONE_MEPS-M; EPS_PA=EPS+A; EPS_PB=EPS+B PI_EPS=M_PI*EPS; EPS_PA_PM=EPS_PA+M; EPS_PB_PM=EPS_PB+M GAMMA_C=ONE/GAMMA_INV(C) GAMMA_INV_EPS_PA_PM=GAMMA_INV(EPS_PA_PM) GAMMA_INV_EPS_PB_PM=GAMMA_INV(EPS_PB_PM) GAMMA_PROD=GAMMA_C*GAMMA_INV_EPS_PA_PM*GAMMA_INV_EPS_PB_PM GAMMA_INV_ONE_MEPS=GAMMA_INV(ONE_MEPS) IF(M.EQ.0) THEN A_TERM=ZERO ELSE IF(INF_NORM(ONE_MEPS_MM & *(LOG(ONE + ABS(ONE_MEPS_MM))-ONE)).LT.300.0d0) THEN A_TERM=GAMMA_PROD*A_SUM_INIT(M,EPS,GAMMA_INV_ONE_MEPS) ELSE A_TERM=EXP(LOG_GAMMA(C)-LOG_GAMMA(EPS_PA_PM)& -LOG_GAMMA(EPS_PB_PM)+LOG_A_SUM_INIT(M,EPS)) IF((AIMAG(A).EQ.ZERO).AND.(AIMAG(B).EQ.ZERO)& .AND.(AIMAG(C).EQ.ZERO)) THEN A_TERM=REAL(A_TERM,wp) ENDIF ENDIF A_SUM=A_TERM POW_MZP1_M = MZP1**M B_TERM=B_SUM_INIT_PS_ONE(A,B,GAMMA_C,GAMMA_INV_ONE_MEPS, & GAMMA_INV_EPS_PA_PM,GAMMA_INV_EPS_PB_PM,MZP1,M,EPS)*POW_MZP1_M PROD_B=POW_MZP1_M DO N=0,M_M2 RATIO=(A+N)*(B+N)/(N+ONE) A_TERM=A_TERM*MZP1*RATIO/(N+ONE_MEPS_MM) A_SUM=A_SUM+A_TERM PROD_B = PROD_B*RATIO ENDDO IF(M.GT.0) THEN PROD_B = PROD_B*(A+M-ONE)*(B+M-ONE)/DBLE(M) ENDIF B_EXTRA_TERM = PROD_B*GAMMA_PROD*GAMMA_INV_ONE_MEPS; B_SUM=B_TERM B_PREC=eps_wp*INF_NORM(B_TERM) CALL CV_POLY_DER_TAB_CALC(A,B,ONE_MEPS_MM,MZP1,CV_POLY1_DER_TAB) CALL CV_POLY_DER_TAB_CALC(EPS_PB_PM,EPS_PA_PM,EPS_PM_P1,MZP1, & CV_POLY2_DER_TAB) MIN_N=MAX(MIN_N_CALC(CV_POLY1_DER_TAB),MIN_N_CALC(CV_POLY2_DER_TAB)) POSSIBLE_FALSE_CV=.TRUE.; N=0 DO WHILE(POSSIBLE_FALSE_CV.OR.(INF_NORM(B_TERM).GT.B_PREC)) N_PM_P1=N+M_P1; N_P1=N+ONE; A_PM_PN=A_PM+N; B_PM_PN=B_PM+N EPS_PM_P1_PN=EPS_PM_P1+N; N_P1_MEPS=ONE_MEPS+N EPS_PM_PN=EPS_PM+N; EPS_PA_PM_PN=EPS_PA_PM+N EPS_PB_PM_PN=EPS_PB_PM+N PROD1=EPS_PA_PM_PN*EPS_PB_PM_PN PROD2=EPS_PM_P1_PN*N_P1 PROD3=A_PM_PN*B_PM_PN B_TERM = MZP1*(B_TERM*PROD1/PROD2+B_EXTRA_TERM*(PROD3/N_PM_P1 & -A_PM_PN-B_PM_PN-EPS+PROD1/N_P1)/(EPS_PM_P1_PN*N_P1_MEPS)) B_SUM=B_SUM+B_TERM B_EXTRA_TERM=B_EXTRA_TERM*MZP1*PROD3/(N_PM_P1*N_P1_MEPS) IF(POSSIBLE_FALSE_CV.AND.(N.GT.MIN_N)) THEN POSSIBLE_FALSE_CV = & (CV_POLY_DER_CALC(CV_POLY1_DER_TAB,DBLE(N)).GT.ZERO).OR. & (CV_POLY_DER_CALC(CV_POLY2_DER_TAB,DBLE(N)).GT.ZERO) ENDIF N=N+1 ENDDO IF(EPS.EQ.ZERO) THEN HYP_PS_ONE=PHASE*(A_SUM+B_SUM) RETURN ELSE HYP_PS_ONE=PHASE*(A_SUM+B_SUM)*PI_EPS/SIN(PI_EPS) RETURN ENDIF END FUNCTION HYP_PS_ONE FUNCTION HYP_PS_INFINITY(A,B,C,Z) ! Calculation of the 2F1 power series ! converging with the 1/z transformation ! ! The formula for F(z) in the 1/z transformation holds: ! F(z) = (-1)^m (pi.eps)/sin (pi.eps) [A(z) + B(z)] ! for eps not equal to zero, ! F(z) = (-1)^m [A(z) + B(z)] for eps = 0 ! where m = |Re(b-a)], eps = b-a-m, ! A(z) = \sum_{n=0}^{m-1} alpha[n] z^{-n}, ! B(z) = \sum_{n=0}^{+oo} beta[n] z^{-n}, and: ! ! alpha[0] = [Gamma_inv(1-m-eps)/eps] Gamma_inv(c-a) ! x Gamma_inv(a+m+eps) Gamma(c) ! [Gamma_inv(1-m-eps)/eps] is calculated in A_sum_init. ! alpha[0] is calculated with log[Gamma] ! if the previous expression might overflow, ! and its imaginary part removed if a, b and c are real. ! alpha[n+1] = (a+n)(1-c+a+n)/[(n+1)(1-m-eps+n)] alpha[n], ! n in [0:m-2]. ! ! beta[0] is defined in B_sum_init_PS_infinity function comments. ! gamma[0] = Gamma(c) (a)_m (1-c+a)_m z^{-m} Gamma_inv(a+m+eps) ! x Gamma_inv(c-a) Gamma_inv(m+1) Gamma_inv(1-eps) ! ! beta[n+1] = (a+m+n+eps)(1-c+a+m+n+eps)/[(m+n+1+eps)(n+1)] beta[n] ! + [(a+m+n)(1-c+a+m+n)/(m+n+1) - (a+m+n) - (1-c+a+m+n) ! - eps + (a+m+n+eps)(1-c+a+m+n+eps)/(n+1)] ! x gamma[n]/[(n+m+1+eps)(n+1+eps)], n >= 0. ! gamma[n+1] = (a+m+n)(b+m+n)/[(m+n+1)(n+1-eps)] gamma[n], n >= 0. ! ! B(z) converges <=> |z| > 1 ! The test of convergence is |beta[n] z^{-n}|oo < 1E-15 |beta[0]|oo ! for n large enough so that false convergence cannot occur. ! ! Variables ! --------- ! a,b,c,z: a,b,c parameters and z argument of 2F1(a,b,c,z) ! m,phase,m_p1,eps,a_mc_p1,one_meps, ! one_meps_pm,a_pm,a_mc_p1_pm,cma: |Re(b-a)], (-1)^m, m+1, b-a-m, ! 1-c+a, 1-eps, 1-eps-m, a+m, 1-c+a+m, c-a ! eps_pa,eps_pm_p1,eps_pa_mc_p1_pm,Pi_eps,eps_pa_pm,eps_pm,Gamma_c: ! eps+a, eps+m+1, eps+1-c+a+m, pi.eps, eps+a+m, eps+m, Gamma(c) ! Gamma_inv_eps_pa_pm,Gamma_inv_cma,z_inv,pow_mz_ma, ! Gamma_inv_one_meps,Gamma_prod: Gamma_inv(eps+a+m), Gamma_inv(c-a), ! 1/z, (-z)^(-a), Gamma_inv(1-eps), ! Gamma(c) Gamma_inv(c-a) Gamma_inv(eps+a+m) ! A_first_term,A_sum,A_term: alpha[0], A(z), alpha[n] z^{-n} ! pow_z_inv_m,B_first_term,prod_B,ratio: z^{-m}, beta[0], ! (a)_m (1-c+a)_m z^{-m}, (a+n)(1-c+a+n)/(n+1) for n in [0:m-2]. ! B_extra_term,B_term,B_sum,B_prec: ! gamma[n], beta[n] z^{-n}, B(z), 1E-15 |beta[0|oo ! cv_poly1_der_tab,cv_poly2_der_tab: P1'(X) and P2'(X) coefficients ! of the potentials derivatives of P1(X) and P2(X) ! defined in cv_poly_der_tab_calc ! with parameters a1 = a, b1 = 1-c+a, c1 = 1-m-eps, z1 = 1/z ! and a2 = b, b2 = eps+1-c+a+m,c2 = eps+m+1, z2 = 1/z. ! min_n: smallest integer after which false convergence cannot occur. ! It is calculated in min_n_calc with both P1'(X) and P2'(X), ! so one takes the largest integer coming from both calculations. ! possible_false_cv: always true if n < min_n. If n >= min_n, ! it is true if P1'(n) > 0 or P2'(n) > 0. ! If n >= min_n and P1'(n) < 0 and P2'(n) < 0, ! it becomes false and remains as such for the rest of the calculation. ! One can then check if |beta[n] z^n|oo < 1E-15 to truncate the series. ! n,n_pm_p1,n_p1,a_pm_pn,a_mc_p1_pm_pn,eps_pm_p1_pn,n_p1_meps, ! eps_pa_pm_pn,eps_pa_mc_p1_pm_pn,eps_pm_pn: ! index of power series, n+m+1, n+1, a+m+n, 1-c+a+m+n, eps+m+n+1, ! n+1-eps, eps+a+m+n, eps+1-c+a+m+n, eps+m+n, ! prod1,prod2,prod3: (eps+a+m+n)(eps+1-c+a+m+n), ! (eps+m+1+n)(n+1), (a+m+n)(1-c+a+m+n) !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: A,B,C,Z INTEGER(i4) :: N,M,PHASE,M_M2,MIN_N,M_P1 REAL(wp) :: B_PREC,N_P1,N_PM_P1 COMPLEX(wp) :: POW_Z_INV_M COMPLEX(wp) :: HYP_PS_INFINITY,Z_INV,RATIO COMPLEX(wp) :: EPS,A_MC_P1,ONE_MEPS,ONE_MEPS_MM,A_PM,A_MC_P1_PM COMPLEX(wp) :: CMA,EPS_PA,EPS_PM_P1,EPS_PA_MC_P1_PM,PI_EPS COMPLEX(wp) :: EPS_PA_PM,EPS_PM,GAMMA_C,GAMMA_INV_CMA,POW_MZ_MA COMPLEX(wp) :: A_SUM,A_TERM COMPLEX(wp) :: GAMMA_INV_EPS_PA_PM,GAMMA_INV_ONE_MEPS COMPLEX(wp) :: PROD_B,B_EXTRA_TERM,B_TERM,B_SUM,PROD1 COMPLEX(wp) :: A_PM_PN,A_MC_P1_PM_PN,EPS_PM_P1_PN,N_P1_MEPS COMPLEX(wp) :: PROD2,PROD3,GAMMA_PROD COMPLEX(wp) :: EPS_PA_PM_PN,EPS_PA_MC_P1_PM_PN,EPS_PM_PN LOGICAL :: POSSIBLE_FALSE_CV REAL(wp) :: CV_POLY1_DER_TAB(0:3),CV_POLY2_DER_TAB(0:3) M=NINT(REAL(B-A,wp)); M_M2=M-2;M_P1=M+1 IF(MOD(M,2).EQ.0) THEN PHASE=1 ELSE PHASE=-1 ENDIF EPS=B-A-M; A_MC_P1=ONE-C+A; ONE_MEPS=ONE-EPS; ONE_MEPS_MM=ONE_MEPS-M A_PM=A+M; A_MC_P1_PM=A_MC_P1+M; CMA=C-A; EPS_PA=EPS+A EPS_PM=EPS+M; EPS_PM_P1=EPS_PM+ONE; EPS_PA_MC_P1_PM=EPS+A_MC_P1_PM PI_EPS=M_PI*EPS; EPS_PA_PM=EPS_PA+M GAMMA_C=ONE/GAMMA_INV(C); GAMMA_INV_EPS_PA_PM = GAMMA_INV(EPS_PA_PM) GAMMA_INV_ONE_MEPS = GAMMA_INV(ONE_MEPS) GAMMA_INV_CMA=GAMMA_INV(CMA); Z_INV=ONE/Z;POW_MZ_MA=(-Z)**(-A) GAMMA_PROD=GAMMA_C*GAMMA_INV_CMA*GAMMA_INV_EPS_PA_PM IF(M.EQ.0) THEN A_TERM=ZERO ELSE IF(INF_NORM(ONE_MEPS_MM & *(LOG(ONE + ABS(ONE_MEPS_MM))-ONE)).LT.300.0d0) THEN A_TERM=GAMMA_PROD*A_SUM_INIT(M,EPS,GAMMA_INV_ONE_MEPS) ELSE A_TERM=EXP(LOG_GAMMA(C)-LOG_GAMMA(CMA)-LOG_GAMMA(B) & + LOG_A_SUM_INIT(M,EPS)) IF((AIMAG(A).EQ.ZERO).AND.(AIMAG(B).EQ.ZERO).AND. & (AIMAG(C).EQ.ZERO)) THEN A_TERM=REAL(A_TERM,wp) ENDIF ENDIF A_SUM=A_TERM POW_Z_INV_M=Z_INV**M B_TERM=B_SUM_INIT_PS_INFINITY(A,C,GAMMA_C,GAMMA_INV_CMA, & GAMMA_INV_ONE_MEPS,GAMMA_INV_EPS_PA_PM,Z,M,EPS)*POW_Z_INV_M PROD_B=POW_Z_INV_M DO N=0,M_M2 RATIO=(A+N)*(A_MC_P1+N)/(N+ONE) A_TERM = A_TERM*Z_INV*RATIO/(N+ONE_MEPS_MM) A_SUM = A_SUM+A_TERM PROD_B = PROD_B*RATIO ENDDO IF (M.GT.0) THEN PROD_B=PROD_B*(A+M-ONE)*(A_MC_P1+M-ONE)/DBLE(M) ENDIF B_EXTRA_TERM = PROD_B*GAMMA_PROD*GAMMA_INV_ONE_MEPS B_SUM=B_TERM B_PREC=eps_wp*INF_NORM(B_TERM) CALL CV_POLY_DER_TAB_CALC(A,A_MC_P1,ONE_MEPS_MM,Z_INV, & CV_POLY1_DER_TAB) CALL CV_POLY_DER_TAB_CALC(B,EPS_PA_MC_P1_PM,EPS_PM_P1, & Z_INV,CV_POLY2_DER_TAB) MIN_N=MAX(MIN_N_CALC(CV_POLY1_DER_TAB),MIN_N_CALC(CV_POLY2_DER_TAB)) POSSIBLE_FALSE_CV=.TRUE.; N=0 DO WHILE(POSSIBLE_FALSE_CV.OR.(INF_NORM(B_TERM).GT.B_PREC)) N_PM_P1=N+M_P1; N_P1=N+ONE; A_PM_PN=A_PM+N A_MC_P1_PM_PN=A_MC_P1_PM+N; EPS_PM_P1_PN=EPS_PM_P1+N N_P1_MEPS=N_P1-EPS; EPS_PA_PM_PN=EPS_PA_PM+N EPS_PA_MC_P1_PM_PN=EPS_PA_MC_P1_PM+N; EPS_PM_PN=EPS_PM+N PROD1=EPS_PA_PM_PN*EPS_PA_MC_P1_PM_PN; PROD2=EPS_PM_P1_PN*N_P1 PROD3=A_PM_PN*A_MC_P1_PM_PN B_TERM = Z_INV*(B_TERM*PROD1/PROD2+B_EXTRA_TERM*(PROD3/N_PM_P1 & -A_PM_PN-A_MC_P1_PM_PN-EPS+PROD1/N_P1) & /(EPS_PM_P1_PN*N_P1_MEPS)) B_SUM=B_SUM+B_TERM B_EXTRA_TERM=B_EXTRA_TERM*Z_INV*PROD3/(N_PM_P1*N_P1_MEPS) IF(POSSIBLE_FALSE_CV.AND.(N.GT.MIN_N)) THEN POSSIBLE_FALSE_CV = (CV_POLY_DER_CALC( & CV_POLY1_DER_TAB,DBLE(N)).GT.ZERO).OR.(& CV_POLY_DER_CALC(CV_POLY2_DER_TAB,DBLE(N)).GT.ZERO) ENDIF N=N+1 ENDDO IF(EPS.EQ.ZERO) THEN HYP_PS_INFINITY=PHASE*POW_MZ_MA*(A_SUM+B_SUM) RETURN ELSE HYP_PS_INFINITY=PHASE*POW_MZ_MA*(A_SUM+B_SUM)*PI_EPS & /SIN(PI_EPS) RETURN ENDIF END FUNCTION HYP_PS_INFINITY FUNCTION HYP_PS_COMPLEX_PLANE_REST(A,B,C,Z) ! Calculation of F(z) in transformation theory missing zones ! of the complex plane with a Taylor series ! ! If z is close to exp(+/- i.pi/3), no transformation in 1-z, z, ! z/(z-1) or combination of them can transform z in a complex number ! of modulus smaller than a given Rmax < 1 . ! Rmax is a radius for which one considers power series summation ! for |z| > Rmax is too slow to be processed. One takes Rmax = 0.9 . ! Nevertheless, for Rmax = 0.9, ! these zones are small enough to be handled ! with a Taylor series expansion around a point z0 close to z ! where transformation theory can be used to calculate F(z). ! One then chooses z0 to be 0.9 z/|z| if |z| < 1, and 1.1 z/|z| ! if |z| > 1, ! so that hyp_PS_zero or hyp_PS_infinity can be used ! (see comments of these functions above). ! For this z0, F(z) = \sum_{n=0}^{+oo} q[n] (z-z0)^n, with: ! q[0] = F(z0), q[1] = F'(z0) = (a b/c) 2F1(a+1,b+1,c+1,z0) ! q[n+2] = [q[n+1] (n (2 z0 - 1) - c + (a+b+c+1) z0) ! + q[n] (a+n)(b+n)/(n+1)]/(z0(1-z0)(n+2)) ! As |z-z0| < 0.1, it converges with around 15 terms, ! so that no instability can occur for moderate a, b and c. ! Convergence is tested ! with |q[n] (z-z0)^n|oo + |q[n+1] (z-z0)^{n+1}|oo. ! Series is truncated when this test is smaller ! than 1E-15 (|q[0]|oo + |q[1] (z-z0)|oo). ! No false convergence can happen here ! as q[n] behaves smoothly for n -> +oo. ! ! Variables ! --------- ! a,b,c,z: a,b,c parameters and z argument of 2F1(a,b,c,z) ! abs_z,is_abs_z_small: |z|, true if |z| < 1 and false if not. ! z0,zc_z0_ratio,z0_term1,z0_term2: 0.9 z/|z| if |z| < 1, ! and 1.1 z/|z| if |z| > 1, (z-z0)/(z0 (1-z0)), ! 2 z0 - 1, c - (a+b+c+1) z0 ! hyp_PS_z0,dhyp_PS_z0,prec: F(z0), F'(z0) calculated with 2F1 ! as F'(z0) = (a b/c) 2F1(a+1,b+1,c+1,z0), ! precision demanded for series truncation ! equal to 1E-15 (|q[0]|oo + |q[1] (z-z0)|oo). ! n,an,anp1,anp2,sum: index of the series, q[n] (z-z0)^n, ! q[n+1] (z-z0)^{n+1}, q[n+2] (z-z0)^{n+2}, ! truncated sum of the power series. !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: A,B,C,Z INTEGER(i4) :: N REAL(wp) :: ABS_Z,PREC COMPLEX(wp) :: HYP_PS_COMPLEX_PLANE_REST COMPLEX(wp) :: Z0,ZC,ZC_Z0_RATIO,Z0_TERM1,Z0_TERM2 COMPLEX(wp) :: HYP_PS_Z0,DHYP_PS_Z0,AN,ANP1,ANP2 ABS_Z=ABS(Z) IF(ABS_Z.LT.ONE) THEN Z0=0.9_wp*Z/ABS_Z; ZC=Z-Z0; ZC_Z0_RATIO=ZC/(Z0*(ONE-Z0)) Z0_TERM1=TWO*Z0 - ONE; Z0_TERM2=C-(A+B+ONE)*Z0 HYP_PS_Z0=HYP_PS_ZERO(A,B,C,Z0) DHYP_PS_Z0=HYP_PS_ZERO(A+ONE,B+ONE,C+ONE,Z0)*A*B/C ELSE Z0=1.1_wp*Z/ABS_Z; ZC=Z-Z0; ZC_Z0_RATIO=ZC/(Z0*(ONE-Z0)) Z0_TERM1=TWO*Z0 - ONE; Z0_TERM2=C-(A+B+ONE)*Z0 HYP_PS_Z0=HYP_PS_INFINITY(A,B,C,Z0) DHYP_PS_Z0=HYP_PS_INFINITY(A+ONE,B+ONE,C+ONE,Z0)*A*B/C ENDIF AN=HYP_PS_Z0;ANP1=ZC*DHYP_PS_Z0;HYP_PS_COMPLEX_PLANE_REST=AN+ANP1 PREC=eps_wp*(INF_NORM(AN)+INF_NORM(ANP1)); N=0 DO WHILE(INF_NORM(AN).GT.PREC) ANP2=ZC_Z0_RATIO*(ANP1*(N*Z0_TERM1-Z0_TERM2)+AN*ZC*(A+N)*(B+N) & /(N+ONE))/(N+TWO) HYP_PS_COMPLEX_PLANE_REST = HYP_PS_COMPLEX_PLANE_REST + ANP2 N=N+1 AN=ANP1 ANP1=ANP2 ENDDO RETURN END FUNCTION HYP_PS_COMPLEX_PLANE_REST RECURSIVE FUNCTION HYP_2F1(A,B,C,Z) RESULT(RES) !! CPC Michel Gauss hypergeometric function \({}_2F_1(a, b; c; z)\). ! ! Calculation of F(z) for arbitrary z using previous routines ! ! Firstly, it is checked if a,b and c are negative integers. ! If neither a nor b is negative integer but c is, ! F(z) is undefined so that the program stops with an error message. ! If a and c are negative integers with c < a, ! or b and c are negative integers with b < a, ! or c is not negative integer integer but a or b is, ! one is in the polynomial case. ! In this case, if |z| < |z/(z-1)| or z = 1, ! hyp_PS_zero is used directly, as then |z| <= 2 ! and no instability arises with hyp_PS_zero ! as long the degree of the polynomial is small (<= 10 typically). ! If not, one uses the transformation ! F(z) = (1-z)^{-a} 2F1(a,c-b,c,z/(z-1)) if a is negative integer ! or F(z) = (1-z)^{-b} 2F1(b,c-a,c,z/(z-1)) if b is negative integer ! along with hyp_PS_zero. ! Indeed, 2F1(a,c-b,c,X) is a polynomial if a is negative integer, ! and so is 2F1(b,c-a,c,X) if b is negative integer, ! so that one has here |z/(z-1)| <= 2 ! and the stability of the method is the same ! as for the |z| < |z/(z-1)| case. ! If one is in the non-polynomial case, one checks if z >= 1. ! If it is, one is the cut of F(z) ! so that z is replaced by z - 10^{-307}i. ! Then, using F(z) = 2F1(b,a,c,z) ! and F(z) = (1-z)^{c-a-b} 2F1(c-a,c-b,c,z), ! one replaces a,b,c parameters by combinations of them ! so that Re(b-a) >= 0 and Re(c-a-b) >= 0. ! Exchanging a and b does not change convergence properties, ! while having Re(c-a-b) >= 0 accelerates it ! (In hyp_PS_zero, t[n] z^n ~ z^n/(n^{c-a-b}) for n -> +oo). ! If |1-z| < 1E-5, one uses hyp_PS_one ! as the vicinity of the singular point z = 1 is treated properly. ! After that, one compares |z| and |z/(z-1)| ! to R in {0.5,0.6,0.7,0.8,0.9}. ! If one of them is smaller than R, ! one uses hyp_PS_zero without transformation ! or with the transformation F(z) = (1-z)^{-a} 2F1(a,c-b,c,z/(z-1)). ! Then, if both of them are larger than 0.9, ! one compares |1/z|, |(z-1)/z|, |1-z| and |1/(1-z)| ! to R in {0.5,0.6,0.7,0.8,0.9}. ! If one of them is found smaller than R, ! with the condition that |c-b|oo < 5 for (z-1)/z transformation, ! |a,b,c|oo < 5 for |1-z| transformation ! and |a,c-b,c|oo < 5 for |1/(1-z)| transformation, ! the corresponding transformation is used. ! If none of them was smaller than 0.9, ! one is in the missing zones of transformation theory ! so that the Taylor series of hyp_PS_complex_plane_rest is used. ! ! Variables ! --------- ! a,b,c,z: a,b,c parameters and z argument of 2F1(a,b,c,z) ! Re_a,Re_b,Re_c,na,nb,nc,is_a_neg_int,is_b_neg_int,is_c_neg_int: ! real parts of a,b,c, closest integers to a,b,c, ! true if a,b,c is negative integers and false if not. ! zm1,z_over_zm1,z_shift: z-1, z/(z-1), z - 10^{-307}i in case z >= 1. ! ab_condition, cab_condition: true if Re(b-a) >= 0 and false if not, ! true if Re(c-a-b) >= 0 and false if not. ! abs_zm1,abz_z,abs_z_inv,abs_z_over_zm1,abs_zm1_inv,abs_zm1_over_z: ! |z-1|, |z|, |1/z|, |z/(z-1)|, |1/(z-1)|, |(z-1)/z| ! are_ac_small: true if |a|oo < 5 and |c|oo < 5, false if not. ! is_cmb_small: true if |c-b|oo < 5, false if not. ! are_abc_small: true if |a|oo < 5, |b|oo < 5 and |c|oo < 5, ! false if not. ! are_a_cmb_c_small: true if |a|oo < 5, |c-b|oo < 5 and |c|oo < 5, ! false if not. ! R_tab,R: table of radii {0.5,0.6,0.7,0.8,0.9}, one of these radii. ! res: returned result !---------------------------------------------------------------------- COMPLEX(wp),INTENT(IN) :: A,B,C,Z INTEGER(i4) :: NA,NB,NC,I REAL(wp) :: RE_A,RE_B,RE_C,ABS_Z,ABS_ZM1,ABS_Z_OVER_ZM1 REAL(wp) :: ABS_ZM1_OVER_Z,ABS_ZM1_INV,R_TABLE(1:5),R,ABS_Z_INV COMPLEX(wp) :: RES,Z_SHIFT COMPLEX(wp) :: Z_OVER_ZM1,ZM1 LOGICAL :: IS_A_NEG_INT,IS_B_NEG_INT,IS_C_NEG_INT LOGICAL :: AB_CONDITION,CAB_CONDITION,ARE_A_CMB_C_SMALL LOGICAL :: IS_CMB_SMALL,ARE_AC_SMALL,ARE_ABC_SMALL ! RE_A=REAL(A,wp); RE_B=REAL(B,wp); RE_C=REAL(C,wp); NA=NINT(RE_A); NB=NINT(RE_B); NC=NINT(RE_C); IS_A_NEG_INT=A.EQ.NA.AND.NA.LE.0 IS_B_NEG_INT=B.EQ.NB.AND.NB.LE.0 IS_C_NEG_INT=C.EQ.NC.AND.NC.LE.0 ZM1=Z-ONE IF(IS_C_NEG_INT) THEN ABS_Z=ABS(Z); Z_OVER_ZM1 = Z/ZM1 ABS_Z_OVER_ZM1=ABS(Z_OVER_ZM1) IF(IS_A_NEG_INT.AND.(NC.LT.NA)) THEN IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN RES=HYP_PS_ZERO(A,B,C,Z) RETURN ELSE RES=((-ZM1)**(-A))*HYP_PS_ZERO(A,C-B,C,Z_OVER_ZM1) RETURN ENDIF ELSE IF(IS_B_NEG_INT.AND.(NC.LT.NB)) THEN IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN RES=HYP_PS_ZERO(A,B,C,Z) RETURN ELSE RES=((-ZM1)**(-B))*HYP_PS_ZERO(B,C-A,C,Z_OVER_ZM1) RETURN ENDIF ELSE STOP '2F1 UNDEFINED' ENDIF ENDIF IF(IS_A_NEG_INT) THEN ABS_Z=ABS(Z); Z_OVER_ZM1 = Z/ZM1 ABS_Z_OVER_ZM1=ABS(Z_OVER_ZM1) IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN RES=HYP_PS_ZERO(A,B,C,Z) RETURN ELSE RES=((-ZM1)**(-A))*HYP_PS_ZERO(A,C-B,C,Z_OVER_ZM1) RETURN ENDIF ELSE IF(IS_B_NEG_INT) THEN ABS_Z=ABS(Z); Z_OVER_ZM1 = Z/ZM1 ABS_Z_OVER_ZM1=ABS(Z_OVER_ZM1) IF((Z.EQ.ONE).OR.(ABS_Z.LT.ABS_Z_OVER_ZM1)) THEN RES=HYP_PS_ZERO(A,B,C,Z) RETURN ELSE RES=((-ZM1)**(-B))*HYP_PS_ZERO(B,C-A,C,Z_OVER_ZM1) RETURN ENDIF ENDIF IF((REAL(Z,wp).GE.ONE).AND.(AIMAG(Z).EQ.ZERO)) THEN Z_SHIFT=CMPLX(REAL(Z,wp),-1.0e-307_wp,wp) RES=HYP_2F1(A,B,C,Z_SHIFT) RETURN ENDIF AB_CONDITION = (RE_B.GE.RE_A) CAB_CONDITION = (RE_C.GE.RE_A + RE_B) IF ((.NOT.AB_CONDITION).OR.(.NOT.CAB_CONDITION)) THEN IF ((.NOT.AB_CONDITION).AND.(CAB_CONDITION)) THEN RES=HYP_2F1(B,A,C,Z) RETURN ELSE IF((.NOT.CAB_CONDITION).AND.(AB_CONDITION)) THEN RES=((-ZM1)**(C-A-B))*HYP_2F1(C-B,C-A,C,Z) RETURN ELSE RES=((-ZM1)**(C-A-B))*HYP_2F1(C-A,C-B,C,Z) RETURN ENDIF ENDIF ABS_ZM1=ABS(ZM1) IF(ABS_ZM1.LT.1.0e-5_wp) THEN RES=HYP_PS_ONE (A,B,C,-ZM1) RETURN ENDIF ABS_Z=ABS(Z); ABS_Z_OVER_ZM1=ABS_Z/ABS_ZM1; ABS_Z_INV=ONE/ABS_Z ABS_ZM1_OVER_Z=ONE/ABS_Z_OVER_ZM1; ABS_ZM1_INV=ONE/ABS_ZM1 IS_CMB_SMALL = INF_NORM(C-B).LT.5.0_wp; ARE_AC_SMALL = (INF_NORM(A).LT.5.0_wp).AND.(INF_NORM(C).LT.5.0_wp) ARE_ABC_SMALL = ARE_AC_SMALL.AND.(INF_NORM(B).LT.5.0_wp) ARE_A_CMB_C_SMALL = ARE_AC_SMALL.AND.IS_CMB_SMALL R_TABLE=[0.5_wp,0.6_wp,0.7_wp,0.8_wp,0.9_wp] DO I=1,5 R=R_TABLE(I) IF(ABS_Z.LE.R) THEN RES=HYP_PS_ZERO (A,B,C,Z) RETURN ENDIF IF(IS_CMB_SMALL.AND.(ABS_Z_OVER_ZM1.LE.R)) THEN RES=((-ZM1)**(-A))*HYP_PS_ZERO (A,C-B,C,Z/ZM1) RETURN ENDIF ENDDO DO I=1,5 R=R_TABLE(I) IF(ABS_Z_INV.LE.R) THEN RES=HYP_PS_INFINITY (A,B,C,Z) RETURN ENDIF IF(IS_CMB_SMALL.AND.(ABS_ZM1_OVER_Z.LE.R)) THEN RES=((-ZM1)**(-A))*HYP_PS_INFINITY (A,C-B,C,Z/ZM1) RETURN ENDIF IF(ARE_ABC_SMALL.AND.(ABS_ZM1.LE.R)) THEN RES=HYP_PS_ONE (A,B,C,-ZM1) RETURN ENDIF IF(ARE_A_CMB_C_SMALL.AND.(ABS_ZM1_INV.LE.R)) THEN RES=((-ZM1)**(-A))*HYP_PS_ONE (A,C-B,C,-ONE/ZM1) RETURN ENDIF ENDDO RES=HYP_PS_COMPLEX_PLANE_REST (A,B,C,Z) RETURN END FUNCTION HYP_2F1 FUNCTION TEST_2F1(A,B,C,Z,F) ! Test of 2F1 numerical accuracy ! using hypergeometric differential equation ! ! If z = 0, F(z) = 1 so that this value is trivially tested. ! To test otherwise if the value of F(z) is accurate, ! one uses the fact that ! z(z-1) F''(z) + (c - (a+b+1) z) F'(z) - a b F(z) = 0. ! If z is not equal to one, a relative precision test is provided ! by |F''(z) + [(c - (a+b+1) z) F'(z) - a b F(z)]/[z(z-1)]|oo ! /(|F(z)|oo + F'(z)|oo + |F''(z)|oo). ! If z is equal to one, one uses |(c - (a+b+1)) F'(z) - a b F(z)|oo ! /(|F(z)|oo + F'(z)|oo + 1E-307). ! F'(z) and F''(z) are calculated using equalities ! F'(z) = (a b/c) 2F1(a+1,b+1,c+1,z) ! and F'(z) = ((a+1)(b+1)/(c+1)) (a b/c) 2F1(a+2,b+2,c+2,z). ! ! Variables ! --------- ! a,b,c,z: a,b,c parameters and z argument of 2F1(a,b,c,z) ! F,dF,d2F: F(z), F'(z) and F''(z) calculated with hyp_2F1 ! using F'(z) = (a b/c) 2F1(a+1,b+1,c+1,z) ! and F'(z) = ((a+1)(b+1)/(c+1)) (a b/c) 2F1(a+2,b+2,c+2,z). !---------------------------------------------------------------------- COMPLEX(wp), INTENT(IN) :: A,B,C,Z REAL(wp) :: TEST_2F1 COMPLEX(wp) :: F,DF,D2F IF(Z.EQ.ZERO) THEN TEST_2F1=INF_NORM(F-ONE) RETURN ELSE IF(Z.EQ.ONE) THEN DF = HYP_2F1(A+ONE,B+ONE,C+ONE,Z)*A*B/C TEST_2F1=INF_NORM((C-(A+B+ONE))*DF-A*B*F) & /(INF_NORM(F)+INF_NORM(DF)+1e-307_wp) RETURN ELSE DF = HYP_2F1(A+ONE,B+ONE,C+ONE,Z)*A*B/C D2F = HYP_2F1(A+TWO,B+TWO,C+TWO,Z)*A*(A+ONE)*B*(B+ONE) & /(C*(C+ONE)) TEST_2F1=INF_NORM(D2F+((C-(A+B+ONE)*Z)*DF-A*B*F)/(Z*(ONE-Z))) & /(INF_NORM(F)+INF_NORM(DF)+INF_NORM(D2F)) RETURN ENDIF END FUNCTION TEST_2F1 end module cpc_michel