From 47964a0c23aff5fec41543ad4398df4d83cfbd96 Mon Sep 17 00:00:00 2001
From: Luca <luca.naterop@bluewin.ch>
Date: Tue, 9 Apr 2019 11:35:50 +0200
Subject: [PATCH] evaluate GPL with all zi = 0

---
 gpl_module.f90 | 22 +++++++++++++++++++++-
 test.f90       | 27 ++++++++++++++++++---------
 2 files changed, 39 insertions(+), 10 deletions(-)

diff --git a/gpl_module.f90 b/gpl_module.f90
index 8bfd057..e722bc5 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 a8e4efc..1d8f5f1 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)
-- 
GitLab