mpl_module.f90 1.93 KB
Newer Older
Luca's avatar
Luca committed
1 2
 
MODULE mpl_module
3 4
  use globals
  use utils
5
  implicit none
Luca's avatar
Luca committed
6

7 8
CONTAINS 

9
  FUNCTION polylog_series(m,x,n_passed) result(res)
10
    ! Computes the classical polylogarithm Li_m(x) using series representation up to order n
11
    integer :: m
12 13
    integer, optional :: n_passed
    complex(kind=prec) :: x, res
14 15
    integer :: i,n
    integer, allocatable :: j(:)
16
    n = merge(n_passed,GPLInfinity,present(n_passed))  
17
    j = (/(i, i=1,n,1)/) 
18 19
    res = sum(x**j / j**m)
  END FUNCTION polylog_series
20

21
  FUNCTION MPL_converges(m,x)
Luca's avatar
Luca committed
22 23
    ! checks if the MPL converges 
    complex(kind=prec) :: x(:)
24
    integer :: m(:)
Luca's avatar
Luca committed
25
    logical :: MPL_converges
26
    MPL_converges = .false.
Luca's avatar
Luca committed
27
    if(abs(product(x)) < 1) then
28 29 30
      if(m(1) /= 1 .or. x(1) /= 1) then
        MPL_converges = .true.
      end if
Luca's avatar
Luca committed
31
    end if
Luca's avatar
Luca committed
32
  END FUNCTION MPL_converges
Luca's avatar
Luca committed
33

Luca's avatar
Luca committed
34
  recursive FUNCTION MPL(m, x, n_passed) result(res)
35 36
    ! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
    integer :: m(:)
Luca's avatar
Luca committed
37 38
    complex(kind=prec) :: x(:)
    complex(kind=prec) :: res
39 40
    integer, optional :: n_passed
    integer :: i, n
41

42
    n = merge(n_passed,GPLInfinity,present(n_passed))  
Luca's avatar
Luca committed
43

44 45 46 47 48 49
    if(size(m) /= size(x)) then
      print*, 'Error: m and x must have the same length'
    end if
    
    if(size(m) == 1) then
      ! base case
50
      res = polylog_series(m(1),x(1), n)
51 52 53 54
    else 
      ! recursion step
      res = 0
      do i = 1, n    
Luca's avatar
Luca committed
55
        res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1)
56 57 58 59
      end do

      ! a nicer way to do it would be but problem is i
      ! i = (/(j, j=1,n, 1)/)
Luca Naterop's avatar
Luca Naterop committed
60 61
      ! res = sum( x(1)**i / i**m(1) * MPL(m(2:), x(2:), i(1) - 1) )
      ! we could get around this problem by rewriting MPL to operate on each i and returning an array
62 63
      
    end if
Luca's avatar
Luca committed
64
  END FUNCTION MPL
Luca's avatar
Luca committed
65

Luca's avatar
Luca committed
66
END MODULE mpl_module
67

Luca's avatar
Luca committed
68
! PROGRAM test
Luca's avatar
Luca committed
69
!   use mpl_module
Luca's avatar
Luca committed
70
!   logical :: result
71
!   result = MPL_converges( dcmplx((/1.0d0,.7d0,.3d0/)), (/ 1,2,1 /) )
Luca's avatar
Luca committed
72 73
!   print*, result
! end PROGRAM test