diff --git a/gpl_module.f90 b/gpl_module.f90 index e1109d605c42a96b890ffcd89b57dcd89935a291..bfb44846f036a955689b9147f3c8509fe6196fb4 100644 --- a/gpl_module.f90 +++ b/gpl_module.f90 @@ -10,6 +10,7 @@ CONTAINS FUNCTION zeta(n) real(kind=prec) :: values(9), zeta integer :: n + if(n == 1) print*, 'WARNING: zeta(1) divergent' values = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, & 1.03692775514337, 1.0173430619844488, 1.008349277381923, & 1.0040773561979441, 1.0020083928260821, 1.000994575127818/) @@ -72,9 +73,11 @@ CONTAINS return end if + ! if depth one and m = 1 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))) + return end if END FUNCTION pending_integral @@ -82,10 +85,12 @@ CONTAINS FUNCTION reduce_to_convergent(a,y2) result(res) complex(kind=prec) :: a(:), y2, res, sr integer :: i + res = 0 i = min_index(abs(a)) sr = a(i) if(i == 1) then + print*, 's_r at at first place' 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]) & @@ -93,25 +98,41 @@ CONTAINS return end if + 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. complex(kind=prec) :: z_flat(:), y, res complex(kind=prec), allocatable :: z(:), s(:,:) - integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i + integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i, m_1 integer, allocatable :: m(:) + logical :: is_depth_one call print_G(z_flat,y) - ! only two arguments? -> we get a logarithm - ! need make convergent? if(.not. is_convergent(z_flat,y)) then + print*, 'need to make convergent' - if(size(z_flat) == 1) then - print*, 'we get simply a log' - res = log(y-z_flat(1)) - log(-z_flat(1)) + + m_prime = get_condensed_m(z_flat) + m_1 = m_prime(1) + is_depth_one = (count((m_prime>0)) == 1) + if(is_depth_one) then + ! case m = 1 + if(size(z_flat) == 1) then + print*, 'case m = 1' + res = log(y-z_flat(1)) - log(-z_flat(1)) + return + end if + + ! case m >= 2 + print*, 'case m > 1' + print*, 'm = ', m_1 + res = -zeta(m_1) & + + pending_integral([y, cmplx(0.0)], m_1-1, [zero_array(m_1-2), y]) & + - pending_integral([z_flat(m_1), cmplx(0.0)], m_1-1, [zero_array(m_1-2), y]) return end if @@ -191,9 +212,9 @@ CONTAINS do i = 1,k x(i) = merge(y/z(1), z(i-1)/z(i),i == 1) end do - print*, 'computed using Li with ' - print*, 'm = ', m - print*, 'x = ', x + ! print*, 'computed using Li with ' + ! print*, 'm = ', m + ! print*, 'x = ', x res = (-1)**k * MPL(m,x) END FUNCTION G_condensed diff --git a/test.f90 b/test.f90 index c9bf88b48f9199fbc9dae7d96f1e87e34a58f63f..625ff8d95bc855a7047c5741f4c4f86c154aaef7 100644 --- a/test.f90 +++ b/test.f90 @@ -13,13 +13,19 @@ PROGRAM TEST complex(kind=prec) :: res real, parameter :: tol = 1.0e-14 logical :: tests_successful = .true. + integer :: i ! call do_MPL_tests() ! call do_GPL_tests() ! call do_shuffle_tests() ! put this somewhere else - res = G_flat(cmplx((/1.0,3.0/)),cmplx(2.0)) - print*, res + ! res = G_flat(cmplx((/0.0,10.0/)),cmplx(20.0)) + ! print*, res + + do i = 1,8 + res = zeta(i) + print*, 'zeta(',i,') =', res + end do ! if(tests_successful) then ! print*, 'All tests passed. ' diff --git a/utils.f90 b/utils.f90 index ac42834716adec4fe70826d648c2b1070b531c3a..3366182aa7a9fc8a154c801189e32c7c18e2fa50 100644 --- a/utils.f90 +++ b/utils.f90 @@ -15,7 +15,7 @@ MODULE utils CONTAINS FUNCTION get_condensed_m(z) result(m) - ! returns condensed m where the ones not needed are filled with 0 + ! returns condensed m where the ones not needed are filled with 0 (returns same size as z) complex(kind=prec), intent(in) :: z(:) integer :: m(size(z)), pos, i m = 1 @@ -87,9 +87,9 @@ CONTAINS real(kind=prec) :: v(:), minimum integer :: min_index, i min_index = 1 - minimum = v(1) + minimum = 1e15 do i = 1,size(v) - if(v(i) < minimum) then + if(v(i) < minimum .and. v(i) > zero) then minimum = v(i) min_index = i end if