diff --git a/gpl_module.f90 b/gpl_module.f90 index 1202bdf4584a4beb9cfad375d2ff0debef2767b5..c215e22cefeac4236248cd408e051c84dea98ea4 100644 --- a/gpl_module.f90 +++ b/gpl_module.f90 @@ -60,48 +60,40 @@ CONTAINS print*, 'G(', abs(z_flat), abs(y), ')' END SUBROUTINE print_G - RECURSIVE FUNCTION pending_integral(p,i,g) result(res) - ! reduces a pending integral given by - ! p = (y1, b1, ..., br) - ! i = position of integration variable within g - ! g = arguments of G-function, without the variable that is integrated over - ! to a G function - complex(kind=prec) :: p(:), g(:), res, y2 - integer :: i - res = 0 - - ! is what we have a G-function? - if(size(g)+1 == i) then - res = G_flat( [p(2:size(p)), g], p(1) ) - return - end if - ! case where we use (64) - if(size(g) == 1 .and. i == 1) then - y2 = g(1) - res = pending_integral(p,2,[sub_ieps(y2)]) - pending_integral(p,2,[cmplx(0.0)]) & - + G_flat(p(2:size(p)),p(1)) * log(-sub_ieps(y2)) - return - end if + RECURSIVE FUNCTION pending_integral(p,i,g) result(res) + ! evaluates a pending integral by reducing it to simpler ones and g functions + complex(kind=prec) :: p(:), g(:), res + integer :: i + + res = 0 + if(i == size(g)+1) then + res = G_flat([p(2:size(p)),g], p(1)) + return + end if + + if(size(g) == 1) then + 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))) + end if END FUNCTION pending_integral - FUNCTION reduce_to_convergent(a, y2) result(res) - complex(kind=prec) :: a(:), y2, res, s_r - integer :: min_i - min_i = min_index(abs(a)) - s_r = a(min_i) + FUNCTION reduce_to_convergent(a,y2) result(res) + complex(kind=prec) :: a(:), y2, res, sr + integer :: i res = 0 - ! case that minimum is at first place - if(min_i == 1) then - res = G_flat([cmplx(0.0), a(2:size(a))], y2) ! first term of (64) - res = res + G_flat([y2], s_r) * G_flat(a(2:size(a)), y2) - res = res + pending_integral( [s_r, a(min_i+1)], min_i, [a(3:size(a)), y2] ) - res = res + G_flat([a(min_i+1)],s_r) * G_flat(a(2:size(a)), y2) + i = min_index(abs(a)) + sr = a(i) + if(i == 1) then + res = G_flat([cmplx(0), a(i+1:size(a))], y2) & + + G_flat([y2], sr) * G_flat(a(i+1:size(a)), y2) & + + pending_integral([sr, a(i+1)], i, [a(i+2:size(a)), y2]) & + - G_flat([a(i+1)], sr) * G_flat(a(i+1:size(a)), y2) return end if - END FUNCTION reduce_to_convergent + END FUNCTION reduce_to_convergent RECURSIVE FUNCTION G_flat(z_flat,y) result(res) ! Calls G function with flat arguments, that is, zeroes not passed through the m's. @@ -154,7 +146,7 @@ CONTAINS res = G_condensed(m,z,y,size(m)) deallocate(m) deallocate(z) - END FUNCTION G_flat + END FUNCTION G_flat RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res) ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) diff --git a/test.f90 b/test.f90 index c4fcacd0e8e3b4cc75cb5ebfb604a24d8419e406..d7198bf5234bb503359aaaa4def969a59ffff4d2 100644 --- a/test.f90 +++ b/test.f90 @@ -18,7 +18,7 @@ PROGRAM TEST ! call do_GPL_tests() ! call do_shuffle_tests() ! put this somewhere else - res = G_flat(cmplx((/0.3,2.2/)),cmplx(2.0)) + res = G_flat(cmplx((/1.0,10.0/)),cmplx(2.0)) print*, res ! if(tests_successful) then