diff --git a/src/maths_functions.f90 b/src/maths_functions.f90 index f3231e3498fdeba12bde8817781733f1c6cdaffa..5f27bde88b2efd74d08715e9b03b65b471f0d119 100644 --- a/src/maths_functions.f90 +++ b/src/maths_functions.f90 @@ -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)