minor bug, polylog zeta function case

src/eval.f90

@@ -12,9 +12,9 @@ PROGRAM eval
 
   call parse_cmd_args()
 
-  res = GPL([0.1,-1.,1.])
+  res = GPL([-1.,0.,0.,0.,1.])
   print*, res
-  
+
   ! res = pending_integral(cmplx([1,3]),2,cmplx([2,4]))
   ! print*, res
 

src/gpl_module.f90

@@ -13,32 +13,6 @@ MODULE gpl_module
 
 CONTAINS 
 
-  FUNCTION zeta(n) 
-    real(kind=prec) :: values(9), zeta
-    integer :: n
-    values = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, &
-               1.03692775514337, 1.0173430619844488, 1.008349277381923, & 
-               1.0040773561979441, 1.0020083928260821, 1.000994575127818/)
-    zeta = values(n-1)
-  END FUNCTION zeta
-
-  FUNCTION GPL_has_convergent_series(m,z,y)
-    ! tests if GPL has a convergent series representation
-    integer :: m(:)
-    complex(kind=prec) :: z(:), y
-    logical :: GPL_has_convergent_series
-
-    GPL_has_convergent_series = .false.
-
-    if(all(abs(y) <= abs(z))) then
-      if(m(1) == 1) then 
-        GPL_has_convergent_series = .true. !(abs((y/z(1) - 1)) < zero)
-      else 
-        GPL_has_convergent_series = .true.
-      end if
-    end if
-  END FUNCTION GPL_has_convergent_series
-
   FUNCTION GPL_zero_zi(l,y)
     ! used to compute the value of GPL when all zi are zero
     integer :: l
@@ -67,7 +41,7 @@ CONTAINS
     if(.not. present(y)) print*, 'G(', abs(z_flat), ')'
   END SUBROUTINE print_G
 
-  FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) result(res)
+  RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) result(res)
     ! here what is passed is not the full a vector, only a1, ..., ak without the trailing zeroes
     integer :: m, i, j, n
     complex(kind=prec) :: a(:), y2, s(m), p(:), res
@@ -127,7 +101,7 @@ CONTAINS
       return
     end if
   
-    a = g(1:size(p)-1)
+    a = g(1:size(g)-1)
     y2 = g(size(g)) 
     m = size(g)  ! actually, size(g)-1+1
 
@@ -288,10 +262,9 @@ CONTAINS
 
     if(verb >= 50) call print_G(z_flat,y)
 
-    ! catch G(1,1)
+
     if(size(z_flat) == 1) then
       if( abs(z_flat(1) - y) <= zero ) then
-        print*, 'catch G(1,1)'
         res = 0
         return
       end if
@@ -429,7 +402,7 @@ CONTAINS
     end  if
 
     ! need make convergent?
-    if(.not. GPL_has_convergent_series(m,z,y)) then
+    if(.not. all(abs(y) <= abs(z))) then
       z_flat = get_flattened_z(m,z)
       res = G_flat(z_flat,y)
       return

src/maths_functions.f90

@@ -12,11 +12,20 @@ CONTAINS
     complex(kind=prec) :: x, res
     integer :: i,n
     integer, allocatable :: j(:)
-    n = 1000
+    n = 30
     j = (/(i, i=1,n,1)/) 
     res = sum(x**j / j**m)
   END FUNCTION naive_polylog
 
+  FUNCTION zeta(n) 
+    real(kind=prec) :: values(9), zeta
+    integer :: n
+    values = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, &
+               1.03692775514337, 1.0173430619844488, 1.008349277381923, & 
+               1.0040773561979441, 1.0020083928260821, 1.000994575127818/)
+    zeta = values(n-1)
+  END FUNCTION zeta
+
   FUNCTION Li2(x)
 
    !! Dilogarithm for arguments x < = 1.0
@@ -255,7 +264,11 @@ CONTAINS
     else if(m == 3) then
       res = trilog(x)
     else
-      res = naive_polylog(m,x)
+      if(abs(x-cmplx(1)) < zero) then
+        res = zeta(m)
+      else 
+        res = naive_polylog(m,x)
+      end if
     end if
   END FUNCTION polylog
 

src/mpl_module.f90

@@ -2,6 +2,7 @@
 MODULE mpl_module
   use globals
   use utils
+  use maths_functions
   implicit none
 
 CONTAINS 

src/test.f90

@@ -130,8 +130,8 @@ CONTAINS
     ref = (-8.403785974849544e-6,0.0)
     call test_one_flat(cmplx([-1./xchen,0.,-1./xchen,1./(xchen*z),1.0]),ref,'3.5')
     
-    ! ref = cmplx((0.4925755847450199,2.6389214054743295))
-    ! call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
+    ref = cmplx((0.4925755847450199,2.6389214054743295))
+    call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
     
     ! ref = cmplx((-0.0015317713178859967,-0.00045003911367000565))
     ! call test_one_flat(cmplx([0.,0.,(1-sqrt(1-z**2))/(xchen*z), 1./(xchen*z),1.]),ref,'3.7')
@@ -177,7 +177,6 @@ CONTAINS
 
   end subroutine do_GPL_tests
   
-  
   ! subroutine do_shuffle_tests() 
   !   complex(kind=prec) :: v(2) = cmplx((/1,2/))
   !   complex(kind=prec) :: w(2) = cmplx((/3,4/))