Commit 47964a0c authored by Luca's avatar Luca

evaluate GPL with all zi = 0

parent bd4c9a82
...@@ -6,16 +6,36 @@ MODULE gpl_module ...@@ -6,16 +6,36 @@ MODULE gpl_module
CONTAINS CONTAINS
RECURSIVE FUNCTION factorial(n) result(res)
integer, intent(in) :: n
integer :: res
res = merge(1,n*factorial(n-1),n==1)
END FUNCTION factorial
FUNCTION GPL_zero_zi(l,y)
! used to compute the value of GPL when all zi are zero
integer :: l
complex(kind=prec) :: y, GPL_zero_zi
GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l
END FUNCTION GPL_zero_zi
FUNCTION GPL(m,z,y,k) FUNCTION GPL(m,z,y,k)
! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y)
! assumes zero arguments expressed through the m's
integer :: m(:), k, i integer :: m(:), k, i
complex(kind=prec) :: z(:), x(k), y, GPL complex(kind=prec) :: z(:), x(k), y, GPL
! first check if we have only zero arguments
if(k == 1 .and. z(1) == 0) then
! for that we assume that only one argument was passed, the rest through m1
GPL = GPL_zero_zi(m(1)-1,y)
return
end if
do i = 1,k do i = 1,k
x(i) = merge(y/z(1), z(i-1)/z(i),i == 1) x(i) = merge(y/z(1), z(i-1)/z(i),i == 1)
end do end do
GPL = (-1)**k * MPL(m,x) GPL = (-1)**k * MPL(m,x)
END FUNCTION GPL END FUNCTION GPL
......
...@@ -11,7 +11,7 @@ PROGRAM TEST ...@@ -11,7 +11,7 @@ PROGRAM TEST
real, parameter :: tol = 1.0e-14 real, parameter :: tol = 1.0e-14
integer :: testnr integer :: testnr
call do_MPL_tests() ! call do_MPL_tests()
call do_GPL_tests() call do_GPL_tests()
CONTAINS CONTAINS
...@@ -70,28 +70,37 @@ CONTAINS ...@@ -70,28 +70,37 @@ CONTAINS
end subroutine check_one_MPL end subroutine check_one_MPL
subroutine do_GPL_tests() subroutine do_GPL_tests()
integer :: m(2), k integer :: m2(2), m1(1), k
complex(kind=prec) :: z(2), y, res, ref complex(kind=prec) :: z2(2), z1(1), y, res, ref
print*, 'doing GPL tests...' print*, 'doing GPL tests...'
testnr = 11 testnr = 11
m = (/ 1,1 /) m2 = (/ 1,1 /)
z = dcmplx((/ 1.3d0, 1.1d0 /)) z2 = dcmplx((/ 1.3d0, 1.1d0 /))
y = 0.4 y = 0.4
k = 2 k = 2
res = GPL(m,z,y,k) res = GPL(m2,z2,y,k)
ref = dcmplx(0.0819393734128676) ref = dcmplx(0.0819393734128676)
call check(res,ref) call check(res,ref)
testnr = 12 testnr = 12
m = (/ 3,2 /) m2 = (/ 3,2 /)
z = dcmplx((/ 1.3d0, 1.1d0 /)) z2 = dcmplx((/ 1.3d0, 1.1d0 /))
y = 0.4 y = 0.4
k = 2 k = 2
res = GPL(m,z,y,k) res = GPL(m2,z2,y,k)
ref = dcmplx(0.01592795952537145) ref = dcmplx(0.01592795952537145)
call check(res,ref) call check(res,ref)
testnr = 13
m1 = (/ 4 /)
z1 = dcmplx((/ 0 /))
y = 1.6
k = 1
res = GPL(m1,z1,y,k)
ref = dcmplx(0.0173042341866201179)
call check(res,ref)
end subroutine do_GPL_tests end subroutine do_GPL_tests
subroutine check_one_GPL(m,z,y,k,ref) subroutine check_one_GPL(m,z,y,k,ref)
......
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