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

8 9
CONTAINS 

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

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

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

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

45 46 47 48 49 50
    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
51
      res = polylog_series(m(1),x(1), n)
52 53 54 55
    else 
      ! recursion step
      res = 0
      do i = 1, n    
Luca's avatar
Luca committed
56
        res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1)
57 58 59 60
      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
61 62
      ! 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
63 64
      
    end if
Luca's avatar
Luca committed
65
  END FUNCTION MPL
Luca's avatar
Luca committed
66

Luca's avatar
Luca committed
67
END MODULE mpl_module
68

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