pending integral bug, add small imag part

src/eval.f90

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

src/globals.f90

@@ -4,7 +4,7 @@ MODULE globals
 
   integer, parameter :: prec = selected_real_kind(15,32)  
   integer, parameter :: GPLInfinity = 30   ! the default outermost expansion order for MPLs
-  real, parameter :: epsilon = 1e-15       ! used for the small imaginary part
+  real, parameter :: epsilon = 1e-20       ! used for the small imaginary part
   real, parameter :: zero = 1e-15          ! values smaller than this count as zero
   real, parameter :: pi = 3.14159265358979323846
 

src/gpl_module.f90

@@ -75,50 +75,74 @@ CONTAINS
     integer :: i, m
     res = 0
 
-    if(verb >= 50) print*, 'evaulating PI with p=',abs(p),'i=',abs(i),'g =',abs(g)   
+    if(verb >= 30) print*, 'evaluating PI with p=',abs(p),'i=',abs(i),'g =',abs(g)   
+
+    y1 = p(1)
+    b = p(2:size(p))
+    a = g(1:size(p)-1)
+    y2 = g(size(g)) 
+    m = size(g)  ! actually, size(g)-1+1
+
+    ! if integration variable is not in G-function
+    if(i == 0 .or. size(g) == 0) then
+      if(verb >= 30) print*, 'only integrals in front, we get G-function'
+      res = G_flat(b,y1)
+      return
+    end if
 
     ! if integration variable at end -> we gat a G function 
     if(i == size(g)+1) then
+      if(verb >= 30) print*, 'is just a G-function'
       res = G_flat([p(2:size(p)),g], p(1))
       return
     end if
   
+
     ! if depth one and m = 1 use my (59)
     if(size(g) == 1) then
-      if(verb >= 30) print*, 'depth one and m = 1'
+      if(verb >= 30) print*, 'case depth one and m = 1'
       res = pending_integral(p,2,[sub_ieps(g(1))]) - pending_integral(p,2,[cmplx(0.0)]) &
         + G_flat(p(2:size(p)), p(1)) * log(-sub_ieps(g(1)))
       return
     end if
   
-    y1 = p(1)
-    b = p(2:size(p))
-    a = g(1:size(p)-1)
-    y2 = g(size(g)) 
     
     ! if depth one and m > 1 use my (60)
-    if(all( abs( g(1:size(g)-1) ) < zero)) then 
-      ! thus depth one and m > 1 
-      
-      if(verb >= 30) print*, 'depth one and m > 1'
-      m = size(g)  ! actually, size(g)-1+1
-      res = -zeta(m)*pending_integral(p,1,zeroes(0)) &
-        + pending_integral([y2,cmplx(0.0)],m-1,[zeroes(m-2),y2])*pending_integral(p,1,zeroes(0)) &
+    if(all( abs( g(1:size(g)-1) ) < zero)) then       
+      if(verb >= 30) print*, 'case depth one and m > 1'
+      if(verb >= 50) then 
+        print*, 'map to'
+        print*, 'zeta(',m,')'
+        print*, 'PI with p=',abs(p),'i=',0,'g =',zeroes(0)
+        print*, 'PI with p=',abs([y2,cmplx(0.0)]),'i=',m-1,'g =',[zeroes(m-2),y2]
+        print*, 'PI with p=',abs(p),'i=',0,'g =',zeroes(0)
+        print*, 'PI with p=',[p,cmplx(0.0)],'i=',m-1,'g =',zeroes(0)
+      end if
+      res = -zeta(m)*pending_integral(p,0,zeroes(0)) &
+        + pending_integral([y2,cmplx(0.0)],m-1,[zeroes(m-2),y2])*pending_integral(p,0,zeroes(0)) &
         - pending_integral([p,cmplx(0.0)],m-1,[zeroes(m-2),y2])
-
       return
     end if
   
-  
-    ! case higher depth, s_r at beginning, use my (68)
+    ! case of higher depth, s_r at beginning, use my (68)
     if(i == 1) then
-      if(verb >= 30) print*, 'higher depth, sr at beginning'
+      if(verb >= 30) print*, 'case higher depth, sr at beginning'
       
-      res = pending_integral(p,1,zeroes(0)) * G_flat([cmplx(0.0),a],y2) & 
-        + pending_integral([p,y2],1,zeroes(0)) * G_flat(a,y2) &
-        + pending_integral([p,a(1)],1,[a,y2]) &
-        - pending_integral([p,a(1)],1,zeroes(0)) * G_flat(a,y2)
-    
+      if(verb >= 50) then
+        print*, 'map to (using 68)'
+        print*, 'PI with p=',abs(p),'i=',0,'g =',zeroes(0) 
+        call print_G([cmplx(0.0),a],y2)
+        print*, 'PI with p=',abs([p,y2]),'i=',0,'g =',zeroes(0) 
+        call print_G(a,y2)
+        print*, 'PI with p=',abs([p,a(1)]),'i=',1,'g =',abs([a(2:size(a)),y2])
+        print*, 'PI with p=',abs([p,a(1)]),'i=',0,'g =',zeroes(0)
+        call print_G(a,y2)
+      end if
+
+      res = pending_integral(p,0,zeroes(0)) * G_flat([cmplx(0.0),a],y2) & 
+        + pending_integral([p,y2],0,zeroes(0)) * G_flat(a,y2) &
+        + pending_integral([p,a(1)],1,[a(2:size(a)),y2]) &
+        - pending_integral([p,a(1)],0,zeroes(0)) * G_flat(a,y2)
         return
     end if
     
@@ -224,6 +248,11 @@ CONTAINS
 
     if(verb >= 50) call print_G(z_flat,y)
 
+    ! add small imaginary part if not there
+    do i = 1,size(z_flat)
+      if(abs(aimag(z_flat(i))) < 1e-25) z_flat(i) = add_ieps(z_flat(i))
+      if(abs(aimag(y)) < 1e-25) y = add_ieps(y)
+    end do
 
     ! is just a logarithm? 
     if(all(abs(z_flat) < zero)) then