mpl_module.f90 2.17 KB
Newer Older
Luca's avatar
Luca committed
1 2
 
MODULE mpl_module
3
  implicit none
Luca's avatar
Luca committed
4 5

  integer, parameter :: prec = selected_real_kind(15,32)  
6
  integer :: GPLInfinity = 30   ! the default n if it is not passed
7 8 9 10 11 12

CONTAINS 

  FUNCTION dilog(x,n)
    ! Computes the dilog Li_2(x) using the series representation up to order n
    integer :: n
Luca's avatar
Luca committed
13
    complex(kind=prec) :: x, dilog
14 15 16 17 18 19 20
    integer :: i
    integer :: j(n)

    j = (/(i, i=1,n,1)/) 
    dilog = sum(x**j / j**2)
  END FUNCTION dilog

Luca's avatar
Luca committed
21
  FUNCTION polylog(m,x,n)
22
    ! Computes the classical polylogarithm Li_m(x) using series representation up to order n
23
    integer :: m
Luca's avatar
Luca committed
24
    complex(kind=prec) :: x, polylog
25 26
    integer :: i,n
    integer, allocatable :: j(:)
27 28 29 30 31

    j = (/(i, i=1,n,1)/) 
    polylog = sum(x**j / j**m)
  END FUNCTION polylog

32
  FUNCTION MPL_converges(m,x)
Luca's avatar
Luca committed
33 34
    ! checks if the MPL converges 
    complex(kind=prec) :: x(:)
35
    integer :: m(:)
Luca's avatar
Luca committed
36
    logical :: MPL_converges
37
    MPL_converges = .false.
Luca's avatar
Luca committed
38
    if(abs(product(x)) < 1) then
39 40 41
      if(m(1) /= 1 .or. x(1) /= 1) then
        MPL_converges = .true.
      end if
Luca's avatar
Luca committed
42
    end if
Luca's avatar
Luca committed
43
  END FUNCTION MPL_converges
Luca's avatar
Luca committed
44

Luca's avatar
Luca committed
45
  recursive FUNCTION MPL(m, x, n_passed) result(res)
46 47
    ! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
    integer :: m(:)
Luca's avatar
Luca committed
48 49
    complex(kind=prec) :: x(:)
    complex(kind=prec) :: res
50 51
    integer, optional :: n_passed
    integer :: i, n
52

53
    n = merge(n_passed,GPLInfinity,present(n_passed))  
Luca's avatar
Luca committed
54

55 56 57 58 59 60 61 62 63 64 65
    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
      res = polylog(m(1),x(1), n)
    else 
      ! recursion step
      res = 0
      do i = 1, n    
Luca's avatar
Luca committed
66
        res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1)
67 68 69 70
      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
71 72
      ! 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
73 74
      
    end if
Luca's avatar
Luca committed
75
  END FUNCTION MPL
Luca's avatar
Luca committed
76

Luca's avatar
Luca committed
77
END MODULE mpl_module
78

Luca's avatar
Luca committed
79
! PROGRAM test
Luca's avatar
Luca committed
80
!   use mpl_module
Luca's avatar
Luca committed
81
!   logical :: result
82
!   result = MPL_converges( dcmplx((/1.0d0,.7d0,.3d0/)), (/ 1,2,1 /) )
Luca's avatar
Luca committed
83 84
!   print*, result
! end PROGRAM test