Commit 2c742a80 authored by Luca's avatar Luca

make order n in multiple polylog and classical polylog optional

parent e2de68e8
MODULE functions
use utils
implicit none
integer :: GPLInfinity = 30 ! the default n if it is not passed
CONTAINS
......@@ -15,28 +17,30 @@ CONTAINS
dilog = sum(x**j / j**2)
END FUNCTION dilog
FUNCTION polylog(m,x,n)
FUNCTION polylog(m,x,n_passed)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m,n
integer :: m
complex :: x, polylog
integer :: i
integer :: j(n)
integer :: i,n
integer, optional :: n_passed
integer, allocatable :: j(:)
n = merge(n_passed,GPLInfinity,present(n_passed))
allocate(j(n))
j = (/(i, i=1,n,1)/)
polylog = sum(x**j / j**m)
deallocate(j)
END FUNCTION polylog
recursive FUNCTION multiple_polylog(m, x, n) result(res)
recursive FUNCTION multiple_polylog(m, x, n_passed) 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
integer, optional :: n_passed
integer :: i, n
n = merge(n_passed,GPLInfinity,present(n_passed))
if(size(m) /= size(x)) then
print*, 'Error: m and x must have the same length'
end if
......@@ -70,15 +74,15 @@ PROGRAM test
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*, dilog((0.8,0),20) ! should be 1.07479 + 0i
! print*, dilog((0.2,0.5),20) ! 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
! print*, polylog(2,(0.2,0.5),20) ! should be 0.133909 + 0.537628i
result = polylog(5,(0.2,0.5),20) ! should be 0.192872345 + 0.505898833i
! result = multiple_polylog((/ 5 /),(/ (0.2,0.5) /), 100)
! result = multiple_polylog((/ 5 /),(/ (0.2,0.5) /),10)
! print*, 'result = ', result
! result = multiple_polylog((/ 5, 5 /),(/ (0.8,0),(0.3,0.5) /), 100)
! result = multiple_polylog((/ 5, 5 /),(/ (0.8,0),(0.3,0.5) /), 20)
! print*, 'result = ', result
END PROGRAM test
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment