Commit d8ec8762 authored by Luca's avatar Luca

remove trailing zeros, call G with flat args

parent a4dc981b
MODULE gpl_module
use mpl_module
use utils
implicit none
CONTAINS
RECURSIVE FUNCTION factorial(n) result(res)
integer, intent(in) :: n
integer :: res
......@@ -43,25 +44,60 @@ CONTAINS
! used to compute the value of GPL when all zi are zero
integer :: l
complex(kind=prec) :: y, GPL_zero_zi
print*, 'computed value using zero'
GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l
END FUNCTION GPL_zero_zi
FUNCTION GPL(m,z,y,k)
FUNCTION G_with_flat_args(z_flat,y) result(res)
complex(kind=prec) :: z_flat(:), y, res
complex(kind=prec), allocatable :: z(:)
integer :: m_prime(size(z_flat)), condensed_size
integer, allocatable :: m(:)
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)
res = GPL(m,z,y,size(m))
deallocate(m)
deallocate(z)
END FUNCTION G_with_flat_args
RECURSIVE FUNCTION GPL(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
complex(kind=prec) :: z(:), x(k), y, GPL
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))
! are all z_i = 0 ?
if(k == 1 .and. z(1) == 0) then
! for that we assume that only one argument was passed, the rest through m1
GPL = GPL_zero_zi(m(1),y)
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)) < zero ) then
print*, 'need to remove trailing zeros'
! flatten z
z_flat = get_flattened_z(m,z)
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
end if
! need make convergent?
if(.not. GPL_has_convergent_series(m,z,y,k)) then
print*, 'need to make convergent'
......@@ -70,12 +106,15 @@ CONTAINS
else
print*, ' ', 'does not have convergent series representation'
end if
res = 0
return
end if
do i = 1,k
x(i) = merge(y/z(1), z(i-1)/z(i),i == 1)
end do
GPL = (-1)**k * MPL(m,x)
print*, 'computed using MPL'
res = (-1)**k * MPL(m,x)
END FUNCTION GPL
......
......@@ -7,7 +7,7 @@ FFLAGS=-fdefault-real-8 -cpp
LD=gfortran
OBJ=mpl_module.o gpl_module.o
OBJ= utils.o mpl_module.o gpl_module.o
# Rules for main fortran files
......
MODULE mpl_module
use utils
implicit none
integer, parameter :: prec = selected_real_kind(15,32)
integer :: GPLInfinity = 30 ! the default n if it is not passed
CONTAINS
......
......@@ -3,6 +3,7 @@
! These tests assume that GPLInfinity = 30
PROGRAM TEST
use utils
use mpl_module
use gpl_module
implicit none
......@@ -11,7 +12,7 @@ PROGRAM TEST
real, parameter :: tol = 1.0e-14
logical :: tests_successful = .true.
call do_MPL_tests()
! call do_MPL_tests()
call do_GPL_tests()
if(tests_successful) then
......@@ -20,7 +21,7 @@ PROGRAM TEST
print*, 'Some tests failed. '
stop 1
end if
CONTAINS
subroutine check(res, ref)
......@@ -83,9 +84,12 @@ CONTAINS
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 = dcmplx(0.0020332632172573974)
call test_one_GPL((/ 1 /),(/(0.7,epsilon)/),(1.1,1.4),1,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 PROGRAM TEST
......@@ -11,10 +11,11 @@
MODULE utils
implicit none
logical :: print_enabled = .true.
logical :: warnings_enabled = .true.
integer, parameter :: prec = selected_real_kind(15,32)
real, parameter :: zero = 1.0e-50
! logical :: print_enabled = .true.
! logical :: warnings_enabled = .true.
CONTAINS
......@@ -37,7 +38,7 @@ CONTAINS
m(pos:) = 0
END FUNCTION get_condensed_m
FUNCTION get_condenced_z(m, z_in) result(z_out)
FUNCTION get_condensed_z(m, z_in) result(z_out)
! returns condensed z vector
integer :: m(:), i, pos
complex(kind=prec) :: z_in(:), z_out(size(m))
......@@ -46,7 +47,7 @@ CONTAINS
pos = pos + m(i)
z_out(i) = z_in(pos)
end do
END FUNCTION get_condenced_z
END FUNCTION get_condensed_z
FUNCTION get_flattened_z(m,z_in) result(z_out)
! returns flattened version of z based on m and z
......@@ -119,9 +120,9 @@ END MODULE utils
! implicit none
! complex(kind=prec), dimension(3) :: a = cmplx((/1,2,3/))
! complex(kind=prec) :: z_flat(7)
! complex(kind=prec) :: z_flat(2)
! complex(kind=prec), allocatable :: z(:)
! integer :: m_prime(7), condensed_size
! integer :: m_prime(2), condensed_size
! integer, allocatable :: m(:)
! complex(kind=prec) :: b(size(a)+1,size(a)+1)
......@@ -130,19 +131,21 @@ END MODULE utils
! ! b = shuffle_with_zero(a)
! ! call print_as_matrix(b)
! ! ! test condensing
! ! z_flat = cmplx((/0,0,0,2,1,0,1/))
! ! m_prime = get_condensed_m(z_flat)
! ! condensed_size = find_first_zero(m_prime)-1
! ! allocate(m(condensed_size))
! ! allocate(z(condensed_size))
! ! m = m_prime(1:condensed_size)
! ! z = get_condenced_z(m,z_flat)
! ! z_flat = get_flattened_z(m,z)
! ! print*, z_flat
! ! print*, m
! ! deallocate(m)
! ! deallocate(z)
! ! 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)
! 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!
Please register or to comment