diff --git a/gpl_module.f90 b/gpl_module.f90 index 8bfd0571d32554aae9727687a3586e97002f20c1..e722bc5664d1b83b1e440060c3452e01faca3ca2 100644 --- a/gpl_module.f90 +++ b/gpl_module.f90 @@ -6,16 +6,36 @@ MODULE gpl_module 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) ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) + ! assumes zero arguments expressed through the m's integer :: m(:), k, i 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 x(i) = merge(y/z(1), z(i-1)/z(i),i == 1) end do - GPL = (-1)**k * MPL(m,x) END FUNCTION GPL diff --git a/test.f90 b/test.f90 index a8e4efc4d28c97de9ebd926348d379acacd0a15c..1d8f5f1817c57e55e9b5beed1c7a312780c968e4 100644 --- a/test.f90 +++ b/test.f90 @@ -11,7 +11,7 @@ PROGRAM TEST real, parameter :: tol = 1.0e-14 integer :: testnr - call do_MPL_tests() + ! call do_MPL_tests() call do_GPL_tests() CONTAINS @@ -70,28 +70,37 @@ CONTAINS end subroutine check_one_MPL subroutine do_GPL_tests() - integer :: m(2), k - complex(kind=prec) :: z(2), y, res, ref + integer :: m2(2), m1(1), k + complex(kind=prec) :: z2(2), z1(1), y, res, ref print*, 'doing GPL tests...' testnr = 11 - m = (/ 1,1 /) - z = dcmplx((/ 1.3d0, 1.1d0 /)) + m2 = (/ 1,1 /) + z2 = dcmplx((/ 1.3d0, 1.1d0 /)) y = 0.4 k = 2 - res = GPL(m,z,y,k) + res = GPL(m2,z2,y,k) ref = dcmplx(0.0819393734128676) call check(res,ref) testnr = 12 - m = (/ 3,2 /) - z = dcmplx((/ 1.3d0, 1.1d0 /)) + m2 = (/ 3,2 /) + z2 = dcmplx((/ 1.3d0, 1.1d0 /)) y = 0.4 k = 2 - res = GPL(m,z,y,k) + res = GPL(m2,z2,y,k) ref = dcmplx(0.01592795952537145) 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 subroutine check_one_GPL(m,z,y,k,ref)