MODULE maths_functions use globals use utils implicit none CONTAINS FUNCTION naive_polylog(m,x) result(res) ! Computes the classical polylogarithm Li_m(x) using series representation up to order n integer :: m complex(kind=prec) :: x, res integer :: i,n integer, allocatable :: j(:) n = 1000 j = (/(i, i=1,n,1)/) res = sum(x**j / j**m) END FUNCTION naive_polylog FUNCTION Li2(x) !! Dilogarithm for arguments x < = 1.0 real (kind=prec):: X,Y,T,S,A,PI3,PI6,ZERO,ONE,HALF,MALF,MONE,MTWO real (kind=prec):: C(0:18),H,ALFA,B0,B1,B2,LI2_OLD real (kind=prec):: Li2 integer :: i DATA ZERO /0.0_prec/, ONE /1.0_prec/ DATA HALF /0.5_prec/, MALF /-0.5_prec/ DATA MONE /-1.0_prec/, MTWO /-2.0_prec/ DATA PI3 /3.289868133696453_prec/, PI6 /1.644934066848226_prec/ DATA C( 0) / 0.4299669356081370_prec/ DATA C( 1) / 0.4097598753307711_prec/ DATA C( 2) /-0.0185884366501460_prec/ DATA C( 3) / 0.0014575108406227_prec/ DATA C( 4) /-0.0001430418444234_prec/ DATA C( 5) / 0.0000158841554188_prec/ DATA C( 6) /-0.0000019078495939_prec/ DATA C( 7) / 0.0000002419518085_prec/ DATA C( 8) /-0.0000000319334127_prec/ DATA C( 9) / 0.0000000043454506_prec/ DATA C(10) /-0.0000000006057848_prec/ DATA C(11) / 0.0000000000861210_prec/ DATA C(12) /-0.0000000000124433_prec/ DATA C(13) / 0.0000000000018226_prec/ DATA C(14) /-0.0000000000002701_prec/ DATA C(15) / 0.0000000000000404_prec/ DATA C(16) /-0.0000000000000061_prec/ DATA C(17) / 0.0000000000000009_prec/ DATA C(18) /-0.0000000000000001_prec/ if(X > 1.00000000001_prec) then print*, 'crashes because Li called with bad arguments' elseif(X > 1.0_prec) then X = 1._prec endif IF(X > 0.999999_prec) THEN LI2_OLD=PI6 Li2 = Real(LI2_OLD,prec) RETURN ELSE IF(abs(x-MONE) < zero) THEN LI2_OLD=MALF*PI6 RETURN END IF T=-X IF(T .LE. MTWO) THEN Y=MONE/(ONE+T) S=ONE A=-PI3+HALF*(LOG(-T)**2-LOG(ONE+ONE/T)**2) ELSE IF(T .LT. MONE) THEN Y=MONE-T S=MONE A=LOG(-T) A=-PI6+A*(A+LOG(ONE+ONE/T)) ELSE IF(T .LE. MALF) THEN Y=(MONE-T)/T S=ONE A=LOG(-T) A=-PI6+A*(MALF*A+LOG(ONE+T)) ELSE IF(T .LT. ZERO) THEN Y=-T/(ONE+T) S=MONE A=HALF*LOG(ONE+T)**2 ELSE IF(T .LE. ONE) THEN Y=T S=ONE A=ZERO ELSE Y=ONE/T S=MONE A=PI6+HALF*LOG(T)**2 END IF H=Y+Y-ONE ALFA=H+H B1=ZERO B2=ZERO DO I = 18,0,-1 B0=C(I)+ALFA*B1-B2 B2=B1 B1=B0 ENDDO LI2_OLD=-(S*(B0-H*B2)+A) ! Artificial conversion Li2 = Real(LI2_OLD,prec) END FUNCTION Li2 RECURSIVE FUNCTION dilog(x) result(res) ! evaluates dilog for any argument complex(kind=prec) :: res complex(kind=prec) :: x if(abs(x) <= 1.0) then if(abs(aimag(x)) < zero ) then res = Li2(real(x)) else res = naive_polylog(2,x) endif else res = -dilog(1/x) - (pi**2) /6 - log(add_ieps(-x))**2 / 2 end if END FUNCTION dilog FUNCTION Li3(x) ! Trilogarithm for arguments x < = 1.0 ! This was hacked from LI2 to also follow C332 ! In theory this could also produce Re[Li [x]] for x>1 real (kind=prec):: X,S,A real (kind=prec):: CA(0:18),HA,ALFAA,BA0,BA1,BA2, YA real (kind=prec):: CB(0:18),HB,ALFAB,BB0,BB1,BB2, YB DATA CA(0) / 0.4617293928601208/ DATA CA(1) / 0.4501739958855029/ DATA CA(2) / -0.010912841952292843/ DATA CA(3) / 0.0005932454712725702/ DATA CA(4) / -0.00004479593219266303/ DATA CA(5) / 4.051545785869334e-6/ DATA CA(6) / -4.1095398602619446e-7/ DATA CA(7) / 4.513178777974119e-8/ DATA CA(8) / -5.254661564861129e-9/ DATA CA(9) / 6.398255691618666e-10/ DATA CA(10) / -8.071938105510391e-11/ DATA CA(11) / 1.0480864927082917e-11/ DATA CA(12) / -1.3936328400075057e-12/ DATA CA(13) / 1.8919788723690422e-13/ DATA CA(14) / -2.6097139622039465e-14/ DATA CA(15) / 3.774985548158685e-15/ DATA CA(16) / -5.671361978114946e-16/ DATA CA(17) / 1.1023848202712794e-16/ DATA CA(18) / -5.0940525990875006e-17/ DATA CB(0) / -0.016016180449195803/ DATA CB(1) / -0.5036424400753012/ DATA CB(2) / -0.016150992430500253/ DATA CB(3) / -0.0012440242104245127/ DATA CB(4) / -0.00013757218124463538/ DATA CB(5) / -0.000018563818526041144/ DATA CB(6) / -2.841735345177361e-6/ DATA CB(7) / -4.7459967908588557e-7/ DATA CB(8) / -8.448038544563037e-8/ DATA CB(9) / -1.5787671270014e-8/ DATA CB(10) / -3.0657620579122164e-9/ DATA CB(11) / -6.140791949281482e-10/ DATA CB(12) / -1.2618831590198e-10/ DATA CB(13) / -2.64931268635803e-11/ DATA CB(14) / -5.664711482422879e-12/ DATA CB(15) / -1.2303909436235178e-12/ DATA CB(16) / -2.7089360852246495e-13/ DATA CB(17) / -6.024075373994343e-14/ DATA CB(18) / -1.2894320641440237e-14/ real (kind=prec):: Li3 real (kind=prec), parameter :: zeta2 = 1.6449340668482264365 real (kind=prec), parameter :: zeta3 = 1.2020569031595942854 integer :: i if(x > 1.00000000001_prec) then print*, 'need to crash Li3, since not convergent' elseif(x > 1.0_prec) then x = 1._prec endif IF(X > 0.999999_prec) THEN LI3=zeta3 RETURN ELSE IF( abs(x+1) < zero) THEN LI3=-0.75_prec*zeta3 RETURN END IF IF(X .LE. -1._prec) THEN YA=1._prec/x ; YB=0._prec S=-1._prec A=-LOG(-X)*(zeta2+LOG(-x)**2/6._prec) ELSE IF(X .LE. 0._prec) THEN YA=x ; YB=0._prec S=-1._prec A=0._prec ELSE IF(X .LE. 0.5_prec) THEN YA=0._prec ; YB=x S=-1._prec A=0._prec ELSE IF(X .LE. 1._prec) THEN YA=(x-1._prec)/x ; YB=1._prec-x S=1._prec A=zeta3 + zeta2*Log(x) - (Log(1._prec - X)*Log(X)**2)/2._prec + Log(X)**3/6._prec ELSE IF(X .LE. 2._prec) THEN YA=1._prec - X ; YB=(X-1._prec)/X S=1._prec A=zeta3 + zeta2*Log(x) - (Log(X - 1._prec)*Log(X)**2)/2._prec + Log(X)**3/6._prec ELSE YA=0._prec ; YB=1._prec/X S=-1._prec A=2*zeta2*Log(x)-Log(x)**3/6._prec END IF HA=-2._prec*YA-1._prec ; HB= 2._prec*YB ALFAA=HA+HA ; ALFAB = HB+HB BA0 = 0. ; BA1=0. ; BA2=0. BB0 = 0. ; BB1=0. ; BB2=0. DO I = 18,0,-1 BA0=CA(I)+ALFAA*BA1-BA2 ; BA2=BA1 ; BA1=BA0 BB0=CB(I)+ALFAB*BB1-BB2 ; BB2=BB1 ; BB1=BB0 ENDDO Li3 = A + S * ( (BA0 - HA*BA2) + (BB0 - HB*BB2) ) END FUNCTION Li3 FUNCTION trilog(x) result(res) ! evaluates trilog for any argument complex(kind=prec) :: res complex(kind=prec) :: x if(abs(x) <= 1.0) then if(abs(aimag(x)) < zero ) then res = Li3(real(x)) else res = naive_polylog(3,x) endif else res = naive_polylog(3,sub_ieps(x)**(-1)) - (log(-sub_ieps(x)))**3/6 - pi**2/6 * log(-sub_ieps(x)) end if END FUNCTION trilog FUNCTION polylog(m,x) result(res) ! computes the polylog integer :: m complex(kind=prec) :: x,res if(verb >= 70) print*, 'called polylog(',m,',',x,')' if(m == 2) then res = dilog(x) else if(m == 3) then res = trilog(x) else res = naive_polylog(m,x) end if END FUNCTION polylog END MODULE maths_functions ! PROGRAM test ! use maths_functions ! implicit none ! complex(kind=prec) :: res ! res = Li3(0.4d0) ! print*, res ! END PROGRAM test