diff --git a/gpl_module.f90 b/gpl_module.f90 index e1eaf4f7e4b899e13f8e6f0c264a0cd715c1bdb0..efecef2e3ac69b80c964b20361887ee9478551a7 100644 --- a/gpl_module.f90 +++ b/gpl_module.f90 @@ -6,7 +6,6 @@ MODULE gpl_module CONTAINS - RECURSIVE FUNCTION factorial(n) result(res) integer, intent(in) :: n integer :: res @@ -37,7 +36,6 @@ CONTAINS GPL_has_convergent_series = .true. end if end if - END FUNCTION GPL_has_convergent_series FUNCTION GPL_zero_zi(l,y) @@ -46,15 +44,39 @@ CONTAINS complex(kind=prec) :: y, GPL_zero_zi print*, 'computed value using zi = 0' GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l - END FUNCTION GPL_zero_zi - FUNCTION G_with_flat_args(z_flat,y) result(res) + 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(:) - integer :: m_prime(size(z_flat)), condensed_size + complex(kind=prec), allocatable :: z(:), s(:,:) + integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i integer, allocatable :: m(:) + + print*, 'G_flat called with args', abs(z_flat) + + ! remove trailing zeroes + k = size(z_flat) + kminusj = find_amount_trailing_zeros(z_flat) + j = k - kminusj + if(all(abs(z_flat) < zero)) then + print*, 'all are zero', abs(z_flat) + res = GPL_zero_zi(k,y) + return + else if(kminusj > 0) then + print*, 'we have',kminusj,'trailing zeroes' + allocate(s(j,j)) + s = shuffle_with_zero(z_flat(1:j-1)) + res = log(y)*G_flat(z_flat(1:size(z_flat)-1),y) + do i = 1,size(s,1) + res = res - G_flat([s(i,:),z_flat(j),zero_array(kminusj-1)], y) + end do + res = res / kminusj + deallocate(s) + return + end if + + ! transform to condensed notation m_prime = get_condensed_m(z_flat) if(find_first_zero(m_prime) == -1) then condensed_size = size(m_prime) @@ -65,40 +87,35 @@ CONTAINS allocate(z(condensed_size)) m = m_prime(1:condensed_size) z = get_condensed_z(m,z_flat) - res = GPL(m,z,y,size(m)) + + res = G_condensed(m,z,y,size(m)) deallocate(m) deallocate(z) - END FUNCTION G_with_flat_args + END FUNCTION G_flat - RECURSIVE FUNCTION GPL(m,z,y,k) result(res) + RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res) ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) ! assumes zero arguments expressed through the m's - integer :: m(:), k, i, kminusj + integer :: m(:), k, i complex(kind=prec) :: z(:), x(k), y, res, c(sum(m)+1,sum(m)+1), z_flat(sum(m)), a(sum(m)-1) - ! print*, 'z = ', abs(get_flattened_z(m,z)) + print*, 'called G_condensed with args' + print*, 'm = ', m + print*, 'z = ', abs(z) + ! are all z_i = 0 ? - if(k == 1 .and. abs(z(1)) == 0) then + if(k == 1 .and. abs(z(1)) < zero) then ! assumes that the zeros at the beginning are passed through m1 res = GPL_zero_zi(m(1),y) return end if - ! need to remove trailing zeros? - if(abs(z(k)) == 0 ) then - print*, 'need to remove trailing zeros' ! which we do in flat form - ! flatten z + ! has trailing zeroes? + if(abs(z(k)) < zero ) then + ! we remove them in flat form z_flat = get_flattened_z(m,z) - res = G_with_flat_args(z_flat,y) - - ! a = z_flat(1: (size(z_flat)-1)) - ! c = shuffle_with_zero(a) - ! res = G_with_flat_args(a,y)*log(y) - ! do i = 2,k - ! res = res - G_with_flat_args(c(i,:),y) - ! end do - ! return + res = G_flat(z_flat,y) end if ! need make convergent? @@ -118,8 +135,7 @@ CONTAINS end do print*, 'computed using MPL' res = (-1)**k * MPL(m,x) - - END FUNCTION GPL + END FUNCTION G_condensed END MODULE gpl_module diff --git a/test.f90 b/test.f90 index fdab131d5163071eaaf26b8c5e61cfe365ab6bbd..a1f9eb4b366b498e64a03c45fc14cd2f4f5910e3 100644 --- a/test.f90 +++ b/test.f90 @@ -67,7 +67,7 @@ CONTAINS character(len=*) :: test_id print*, ' ', 'testing GPL ', test_id, ' ...' - res = GPL(m,z,y,k) + res = G_condensed(m,z,y,k) call check(res,ref) end subroutine test_one_GPL @@ -76,17 +76,18 @@ CONTAINS complex(kind=prec), parameter :: epsilon = 1E-14 print*, 'doing GPL tests...' - ref = dcmplx(0.0819393734128676) - call test_one_GPL((/ 1,1 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.1') + ! ref = dcmplx(0.0819393734128676) + ! call test_one_GPL((/ 1,1 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.1') - ref = dcmplx(0.01592795952537145) - call test_one_GPL((/ 3,2 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.2') + ! ref = dcmplx(0.01592795952537145) + ! call test_one_GPL((/ 3,2 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.2') - ref = dcmplx(0.0020332632172573974) - call test_one_GPL((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3') + ! ref = dcmplx(0.0020332632172573974) + ! call test_one_GPL((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3') ref = dcmplx(0.0020332632172573974) - call test_one_GPL((/ 1,1 /),cmplx((/ 1.7,0.0 /)),cmplx(1.1),2,ref,'2.4') + ref = G_flat(cmplx((/1.7,0.0,0.0/)),cmplx(1.1)) + ! call test_one_GPL((/1,1,1/),cmplx((/ 0.0,1.7,0.0 /)),cmplx(1.1),3,ref,'2.4') end subroutine do_GPL_tests diff --git a/utils.f90 b/utils.f90 index 6d929e73f3416ca0105bf61177e1beaf3e7e3c61..dba4ba2d38361f8013512fcd2e97323f96a0d7f6 100644 --- a/utils.f90 +++ b/utils.f90 @@ -12,7 +12,7 @@ MODULE utils implicit none integer, parameter :: prec = selected_real_kind(15,32) - + real :: zero = 1e-15 ! logical :: print_enabled = .true. ! logical :: warnings_enabled = .true. @@ -20,11 +20,12 @@ CONTAINS FUNCTION get_condensed_m(z) result(m) ! returns condensed m where the ones not needed are filled with 0 - complex(kind=prec) :: z(:), m(size(z)) - integer :: pos = 1, i + complex(kind=prec), intent(in) :: z(:) + integer :: m(size(z)), pos, i m = 1 + pos = 1 do i = 1,size(z) - if(z(i) == 0) then + if(abs(z(i)) < zero) then if(i == size(z)) then pos = pos + 1 else @@ -65,7 +66,7 @@ CONTAINS integer :: res, i res = 0 do i = size(z), 1, -1 - if( z(i) == 0 ) then + if( abs(z(i)) < zero ) then res = res + 1 else exit @@ -91,10 +92,16 @@ CONTAINS integer :: s(2), i s = shape(m) do i = 1,s(1) - print*, m(i,:) + print*, abs(m(i,:)) end do END SUBROUTINE print_as_matrix + FUNCTION zero_array(n) result(res) + integer :: n + complex(kind=prec) :: res(n) + res = 0 + END FUNCTION zero_array + FUNCTION shuffle_with_zero(a) result(res) ! rows of result are shuffles of a with 0 complex :: a(:) @@ -130,35 +137,21 @@ END MODULE utils ! PROGRAM test ! use utils ! implicit none - - -! ! complex(kind=prec), dimension(5) :: a = cmplx((/1,2,3/)) -! ! complex(kind=prec) :: z_flat(2) -! ! complex(kind=prec), allocatable :: z(:) -! ! integer :: m_prime(2), condensed_size -! ! integer, allocatable :: m(:) -! ! complex(kind=prec) :: b(size(a)+1,size(a)+1) - -! ! ! ! test shuffling -! ! ! b = 1 -! ! ! b = shuffle_with_zero(a) -! ! ! call print_as_matrix(b) - -! ! ! test condensing -! ! z_flat = cmplx((/4,0/)) -! ! m_prime = get_condensed_m(z_flat) -! ! if(find_first_zero(m_prime) == -1) then -! ! condensed_size = size(m_prime) -! ! else -! ! condensed_size = find_first_zero(m_prime)-1 -! ! end if -! ! allocate(m(condensed_size)) -! ! allocate(z(condensed_size)) -! ! m = m_prime(1:condensed_size) -! ! z = get_condensed_z(m,z_flat) -! ! z_flat = get_flattened_z(m,z) -! ! deallocate(m) -! ! deallocate(z) + +! complex(kind=prec) :: z(4) +! integer :: m_prime(4), condensed_size +! z = cmplx((/0.0,1.7,0.0,0.0/)) + +! ! transform to condensed notation +! m_prime = get_condensed_m(z) +! print*, abs(z) +! m_prime = get_condensed_m(z) +! print*, abs(z) +! if(find_first_zero(m_prime) == -1) then +! condensed_size = size(m_prime) +! else +! condensed_size = find_first_zero(m_prime)-1 +! end if +! print*, condensed_size ! END PROGRAM test -