Commit e769b3af authored by ulrich_y's avatar ulrich_y

Added expansion for PolyLog around log(z)~0 [Crandall 2006]

parent 94ece98d
......@@ -92,6 +92,58 @@ CONTAINS
END FUNCTION
FUNCTION logz_polylog(n, z) result(res)
! Computes the classical polylogarithm Li_m(z) using series
! representation in log(z). valid for log z < 2pi.
!
! The algorithm works by using (1.4) or [Crandall 2006]
!
! PolyLog[n, z] = Ssum[Zeta[n-m] Log[z]^m / m!, {m,0,Infinity}]
! + Log[z]^(n-1) (HarmonicNumber[n-1] - Log[-Log[z]]) / (n-1)!
!
! where Ssum[..] excludes the singular Zeta[1] term at m = n-1. In
! Fortran, we split the Ssum in a sum from 0..n-2 with positive
! arguments in the Zeta function. The next term m=n we do manually
! to not have to implement Zeta[0] = -1/2 and then we use
! Zeta[n-m] = (-1)**(m-n) * bernoullinumber(1+m-n) / (1+m-n)
! for the remaining terms.
!
! References
! R. E. Crandall, Note on fast polylogarithm computation,
! www.reed.edu/~crandall/papers/Polylog.pdf, January 2006.
! or http://functions.wolfram.com/10.08.06.0024.01
integer, parameter :: nmax = 40
real(kind=prec) :: fac, zetamn
complex(kind=prec), intent(in) :: z
complex(kind=prec) :: res, logz
integer, intent(in) :: n
integer m
logz = log(z)
! The factorial will become a problem later. We do the first few
! terms and iterate later.
fac = real(factorial(n-1),kind=prec)
! Non-sum term
res = logz**(n-1) / fac * (harmonicnumber(n-1) - log(-logz))
! positive arguments in the Zeta function, 0..,n-2
do m=0,n-2
res = res + zeta(n-m) / factorial(m) * logz**m
enddo
! Zeta[0], i.e. m=n case
res = res - 0.5_prec * logz**n / fac / n
! All remaining terms
do m=n+1,n+nmax-1,2
zetamn = (-1)**(m-n) * bernoullinumber(1+m-n) / (1+m-n)
fac = fac * m * (m-1)
res = res + zetamn / fac * logz**m
enddo
END FUNCTION logz_polylog
FUNCTION Li2(x)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment