mpl_module.f90 1.25 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 MPL_converges(m,x)
Luca's avatar
Luca committed
11 12
    ! checks if the MPL converges 
    complex(kind=prec) :: x(:)
13
    integer :: m(:)
Luca's avatar
Luca committed
14
    logical :: MPL_converges
15
    MPL_converges = .false.
Luca's avatar
Luca committed
16
    if(abs(product(x)) < 1) then
17
      if(m(1) /= 1 .or. abs(x(1) - 1) < zero) then
18 19
        MPL_converges = .true.
      end if
Luca's avatar
Luca committed
20
    end if
Luca's avatar
Luca committed
21
  END FUNCTION MPL_converges
Luca's avatar
Luca committed
22

Luca's avatar
Luca committed
23
  recursive FUNCTION MPL(m, x, n_passed) result(res)
24 25
    ! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
    integer :: m(:)
Luca's avatar
Luca committed
26 27
    complex(kind=prec) :: x(:)
    complex(kind=prec) :: res
28 29
    integer, optional :: n_passed
    integer :: i, n
30

31
    n = merge(n_passed,MPLInfinity,present(n_passed))  
Luca's avatar
Luca committed
32

33 34 35
    if(size(m) /= size(x)) then
      print*, 'Error: m and x must have the same length'
    end if
Luca Naterop's avatar
Luca Naterop committed
36 37 38 39

    res = 0


40
    if(size(m) == 1) then
Luca Naterop's avatar
Luca Naterop committed
41 42 43 44 45 46
      ! base case, we get a polylog
      do i = 1, n
    if(i**m(1) < 0) return ! roll over
    if(abs(x(1)**i) < 1.e-250) return
        res = res + x(1)**i / i**m(1)
      end do
47 48 49
    else 
      ! recursion step
      do i = 1, n    
Luca's avatar
Luca committed
50
        res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1)
51 52 53
      end do
      
    end if
Luca's avatar
Luca committed
54
  END FUNCTION MPL
Luca's avatar
Luca committed
55

Luca's avatar
Luca committed
56
END MODULE mpl_module