utils.f90 2.14 KB
Newer Older
Luca's avatar
init  
Luca committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87

! Contains some functions that might be useful later

MODULE utils
  implicit none

CONTAINS
  
  FUNCTION dilog(x,n)
    ! Computes the dilog Li_2(x) using the series representation up to order n
    integer :: n
    complex :: x, dilog
    integer :: i
    integer :: j(n)

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

  FUNCTION polylog(m,x,n)
    ! Computes the classical polylogarithm Li_m(x) using series representation up to order n
    integer :: m,n
    complex :: x, polylog
    integer :: i
    integer :: j(n)
    j = (/(i, i=1,n,1)/) 
    polylog = sum(x**j / j**m)

  END FUNCTION polylog

  recursive FUNCTION multiple_polylog(m, x, n) result(res)
    ! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
    integer :: m(:)
    complex :: x(:)
    complex :: res
    integer :: n, i

    ! print*, 'm = ', m
    ! print*, 'x = ', x
    ! print*, 'n = ', n


    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
      ! print*, 'landed in base case'
      res = polylog(m(1),x(1), n)
    else 
      ! recursion step
      ! print*, 'landed in step'
      res = 0
      do i = 1, n    
        res = res + x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i - 1)
      end do

      ! a nicer way to do it would be but problem is i
      ! i = (/(j, j=1,n, 1)/)
      ! res = sum( x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i(1) - 1) )
    
    end if
  END FUNCTION multiple_polylog

END MODULE utils

PROGRAM test
  ! Used to test the procedures defined in this module
  use utils
  implicit none
  complex :: result

  ! print*, dilog((0.8,0),100)     ! should be 1.07479 + 0i 
  ! print*, dilog((0.2,0.5),100)   ! should be 0.133909 + 0.537628i

  ! print*, polylog(2,(0.2,0.5),100)   ! should be 0.133909 + 0.537628i 
  ! print*, polylog(5,(0.2,0.5),100)   ! should be 0.192872345 + 0.505898833i

  result = multiple_polylog((/ 5 /),(/ (0.2,0.5) /), 100)
  print*, 'result = ', result
  result = multiple_polylog((/ 5, 5 /),(/ (0.8,0),(0.3,0.5) /), 100)
  print*, 'result = ', result

END PROGRAM test