Commit 1be3dd9b by Luca Naterop

### function to count amount of trailing zeroes

parent d8ec8762
 ... @@ -44,12 +44,13 @@ CONTAINS ... @@ -44,12 +44,13 @@ CONTAINS ! used to compute the value of GPL when all zi are zero ! used to compute the value of GPL when all zi are zero integer :: l integer :: l complex(kind=prec) :: y, GPL_zero_zi complex(kind=prec) :: y, GPL_zero_zi print*, 'computed value using zero' print*, 'computed value using zi = 0' GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l END FUNCTION GPL_zero_zi END FUNCTION GPL_zero_zi FUNCTION G_with_flat_args(z_flat,y) result(res) FUNCTION G_with_flat_args(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) :: z_flat(:), y, res complex(kind=prec), allocatable :: z(:) complex(kind=prec), allocatable :: z(:) integer :: m_prime(size(z_flat)), condensed_size integer :: m_prime(size(z_flat)), condensed_size ... @@ -73,29 +74,31 @@ CONTAINS ... @@ -73,29 +74,31 @@ CONTAINS ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) ! assumes zero arguments expressed through the m's ! assumes zero arguments expressed through the m's integer :: m(:), k, i integer :: m(:), k, i, kminusj complex(kind=prec) :: z(:), x(k), y, res, c(sum(m)+1,sum(m)+1), z_flat(sum(m)), a(sum(m)-1) 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*, 'z = ', abs(get_flattened_z(m,z)) ! are all z_i = 0 ? ! are all z_i = 0 ? if(k == 1 .and. abs(z(1)) < zero) then if(k == 1 .and. abs(z(1)) == 0) then ! assumes that the zeros at the beginning are passed through m1 ! assumes that the zeros at the beginning are passed through m1 res = GPL_zero_zi(m(1),y) res = GPL_zero_zi(m(1),y) return return end if end if ! need to remove trailing zeros? ! need to remove trailing zeros? if(abs(z(k)) < zero ) then if(abs(z(k)) == 0 ) then print*, 'need to remove trailing zeros' print*, 'need to remove trailing zeros' ! which we do in flat form ! flatten z ! flatten z z_flat = get_flattened_z(m,z) z_flat = get_flattened_z(m,z) a = z_flat(1: (size(z_flat)-1)) res = G_with_flat_args(z_flat,y) c = shuffle_with_zero(a) res = G_with_flat_args(a,y)*log(y) ! a = z_flat(1: (size(z_flat)-1)) do i = 2,k ! c = shuffle_with_zero(a) res = res - G_with_flat_args(c(i,:),y) ! res = G_with_flat_args(a,y)*log(y) end do ! do i = 2,k return ! res = res - G_with_flat_args(c(i,:),y) ! end do ! return end if end if ! need make convergent? ! need make convergent? ... ...
 ... @@ -85,11 +85,9 @@ CONTAINS ... @@ -85,11 +85,9 @@ CONTAINS ref = dcmplx(0.0020332632172573974) ref = dcmplx(0.0020332632172573974) call test_one_GPL((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3') call test_one_GPL((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3') ! ref = dcmplx(0.0020332632172573974) ref = dcmplx(0.0020332632172573974) ! call test_one_GPL((/ 1,1 /),cmplx((/ 1.7,0.0 /)),cmplx(1.1),2,ref,'2.4') call test_one_GPL((/ 1,1 /),cmplx((/ 1.7,0.0 /)),cmplx(1.1),2,ref,'2.4') ! ref = dcmplx(0.0020332632172573974) ! call test_one_GPL((/ 1 /),(/(0.7,epsilon)/),(1.1,1.4),1,ref,'2.4') end subroutine do_GPL_tests end subroutine do_GPL_tests END PROGRAM TEST END PROGRAM TEST
 ... @@ -12,7 +12,6 @@ ... @@ -12,7 +12,6 @@ MODULE utils MODULE utils implicit none implicit none integer, parameter :: prec = selected_real_kind(15,32) integer, parameter :: prec = selected_real_kind(15,32) real, parameter :: zero = 1.0e-50 ! logical :: print_enabled = .true. ! logical :: print_enabled = .true. ! logical :: warnings_enabled = .true. ! logical :: warnings_enabled = .true. ... @@ -61,6 +60,19 @@ CONTAINS ... @@ -61,6 +60,19 @@ CONTAINS end do end do END FUNCTION get_flattened_z END FUNCTION get_flattened_z FUNCTION find_amount_trailing_zeros(z) result(res) complex(kind=prec) :: z(:) integer :: res, i res = 0 do i = size(z), 1, -1 if( z(i) == 0 ) then res = res + 1 else exit end if end do END FUNCTION find_amount_trailing_zeros FUNCTION find_first_zero(v) result(res) FUNCTION find_first_zero(v) result(res) ! returns index of first zero, or -1 if there is no zero ! returns index of first zero, or -1 if there is no zero integer :: v(:), i, res integer :: v(:), i, res ... @@ -119,33 +131,34 @@ END MODULE utils ... @@ -119,33 +131,34 @@ END MODULE utils ! use utils ! use utils ! implicit none ! implicit none ! complex(kind=prec), dimension(3) :: a = cmplx((/1,2,3/)) ! complex(kind=prec) :: z_flat(2) ! ! complex(kind=prec), dimension(5) :: a = cmplx((/1,2,3/)) ! complex(kind=prec), allocatable :: z(:) ! ! complex(kind=prec) :: z_flat(2) ! integer :: m_prime(2), condensed_size ! ! complex(kind=prec), allocatable :: z(:) ! integer, allocatable :: m(:) ! ! integer :: m_prime(2), condensed_size ! complex(kind=prec) :: b(size(a)+1,size(a)+1) ! ! integer, allocatable :: m(:) ! ! complex(kind=prec) :: b(size(a)+1,size(a)+1) ! ! ! test shuffling ! ! b = 1 ! ! ! ! test shuffling ! ! b = shuffle_with_zero(a) ! ! ! b = 1 ! ! call print_as_matrix(b) ! ! ! b = shuffle_with_zero(a) ! ! ! call print_as_matrix(b) ! ! test condensing ! z_flat = cmplx((/4,0/)) ! ! ! test condensing ! m_prime = get_condensed_m(z_flat) ! ! z_flat = cmplx((/4,0/)) ! if(find_first_zero(m_prime) == -1) then ! ! m_prime = get_condensed_m(z_flat) ! condensed_size = size(m_prime) ! ! if(find_first_zero(m_prime) == -1) then ! else ! ! condensed_size = size(m_prime) ! condensed_size = find_first_zero(m_prime)-1 ! ! else ! end if ! ! condensed_size = find_first_zero(m_prime)-1 ! allocate(m(condensed_size)) ! ! end if ! allocate(z(condensed_size)) ! ! allocate(m(condensed_size)) ! m = m_prime(1:condensed_size) ! ! allocate(z(condensed_size)) ! z = get_condensed_z(m,z_flat) ! ! m = m_prime(1:condensed_size) ! z_flat = get_flattened_z(m,z) ! ! z = get_condensed_z(m,z_flat) ! deallocate(m) ! ! z_flat = get_flattened_z(m,z) ! deallocate(z) ! ! deallocate(m) ! ! deallocate(z) ! END PROGRAM test ! END PROGRAM test
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!