Commit 0e52465a authored by Luca's avatar Luca

tests cleanup and convergence criterion of GPL

parent 47964a0c
......@@ -12,6 +12,24 @@ CONTAINS
res = merge(1,n*factorial(n-1),n==1)
END FUNCTION factorial
FUNCTION GPL_has_convergent_series(m,z,y,k)
! tests if GPL has a convergent series representation
integer :: m(:), k
complex(kind=prec) :: z(:), y
logical :: GPL_has_convergent_series
GPL_has_convergent_series = .false.
if(all(abs(y) <= abs(z))) then
if(m(1) == 1) then
GPL_has_convergent_series = (y/z(1) /= 1)
else
GPL_has_convergent_series = .true.
end if
end if
END FUNCTION GPL_has_convergent_series
FUNCTION GPL_zero_zi(l,y)
! used to compute the value of GPL when all zi are zero
integer :: l
......@@ -27,12 +45,18 @@ CONTAINS
integer :: m(:), k, i
complex(kind=prec) :: z(:), x(k), y, GPL
! first check if we have only zero arguments
! 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)-1,y)
return
end if
! do they have convergent series rep?
if(.not. GPL_has_convergent_series(m,z,y,k)) then
print*, ' ', 'does not have convergent series representation'
end if
do i = 1,k
x(i) = merge(y/z(1), z(i-1)/z(i),i == 1)
end do
......
......@@ -29,14 +29,16 @@ CONTAINS
polylog = sum(x**j / j**m)
END FUNCTION polylog
FUNCTION MPL_converges(x)
FUNCTION MPL_converges(m,x)
! checks if the MPL converges
complex(kind=prec) :: x(:)
integer :: m(:)
logical :: MPL_converges
MPL_converges = .false.
if(abs(product(x)) < 1) then
MPL_converges = .true.
else
MPL_converges = .false.
if(m(1) /= 1 .or. x(1) /= 1) then
MPL_converges = .true.
end if
end if
END FUNCTION MPL_converges
......@@ -77,6 +79,6 @@ END MODULE mpl_module
! PROGRAM test
! use mpl_module
! logical :: result
! result = MPL_converges( dcmplx((/10.1d0,.7d0,.3d0/)) )
! result = MPL_converges( dcmplx((/1.0d0,.7d0,.3d0/)), (/ 1,2,1 /) )
! print*, result
! end PROGRAM test
......@@ -9,26 +9,43 @@ PROGRAM TEST
complex(kind=prec) :: res
real, parameter :: tol = 1.0e-14
integer :: testnr
logical :: tests_successful = .true.
! call do_MPL_tests()
call do_MPL_tests()
call do_GPL_tests()
CONTAINS
if(tests_successful) then
print*, 'All tests passed. '
else
print*, 'Some tests failed. '
stop 1
end if
CONTAINS
subroutine check(res, ref)
complex(kind=prec) :: res, ref
real(kind=prec) :: delta
delta = abs((res-ref)/ref)
if(delta < tol) then
print*, testnr,'passed with delta = ', delta
print*, ' ',' passed with delta = ', delta
else
print*, testnr,'not passed with delta = ', delta
stop 1
print*, ' ',' failed with delta = ', delta
tests_successful = .false.
end if
end subroutine check
subroutine test_one_MPL(m,x,ref, test_id)
integer :: m(:)
complex(kind=prec) :: x(:), ref
character(len=*) :: test_id
print*, ' ', 'testing MPL ', test_id, ' ...'
res = MPL(m,x)
call check(res,ref)
end subroutine test_one_MPL
subroutine do_MPL_tests()
integer :: m2(2), m3(3)
integer :: bla
......@@ -36,80 +53,42 @@ CONTAINS
complex(kind=prec) :: res, ref
print*, 'doing MPL tests...'
testnr = 1
m2 = (/ 1,1 /)
x2 = dcmplx((/ 0.3156498673740053, 0.3431255827785649 /))
res = MPL(m2,x2)
ref = dcmplx(0.022696600480693277651633)
call check(res,ref)
call test_one_MPL((/ 1,1 /),cmplx((/ 0.3156498673740053, 0.3431255827785649 /)),ref, '1.1')
testnr = 2
m2 = (/ 1,1 /)
x2 = dcmplx((/ 0.03, 0.5012562893380046 /))
res = MPL(m2,x2)
ref = dcmplx(0.00023134615630308335448329926098409)
call check(res,ref)
call test_one_MPL((/ 1,1 /),cmplx((/ 0.03, 0.5012562893380046 /)),ref, '1.2')
testnr = 3
m3 = (/ 2,1,2 /)
x3 = dcmplx((/ 0.03, 0.5012562893380046, 55.3832 /))
res = MPL(m3,x3)
ref = dcmplx(0.000023446106415452030937059124671151)
call check(res,ref)
call test_one_MPL((/ 2,1,2 /),cmplx((/ 0.03, 0.5012562893380046, 55.3832 /)),ref, '1.3')
end subroutine do_MPL_tests
subroutine check_one_MPL(m,x,ref)
! checks MPL(m,x) against reference value ref
integer :: m(:)
complex(kind=prec) :: x(:)
complex(kind=prec) :: res, ref
subroutine test_one_GPL(m,z,y,k,ref,test_id)
integer :: m(:), k
complex(kind=prec) :: z(:), y, res, ref
character(len=*) :: test_id
res = MPL(m,x)
print*, ' ', 'testing GPL ', test_id, ' ...'
res = GPL(m,z,y,k)
call check(res,ref)
end subroutine check_one_MPL
end subroutine test_one_GPL
subroutine do_GPL_tests()
integer :: m2(2), m1(1), k
complex(kind=prec) :: z2(2), z1(1), y, res, ref
print*, 'doing GPL tests...'
testnr = 11
m2 = (/ 1,1 /)
z2 = dcmplx((/ 1.3d0, 1.1d0 /))
y = 0.4
k = 2
res = GPL(m2,z2,y,k)
ref = dcmplx(0.0819393734128676)
call check(res,ref)
call test_one_GPL((/ 1,1 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.1')
testnr = 12
m2 = (/ 3,2 /)
z2 = dcmplx((/ 1.3d0, 1.1d0 /))
y = 0.4
k = 2
res = GPL(m2,z2,y,k)
ref = dcmplx(0.01592795952537145)
call check(res,ref)
call test_one_GPL((/ 3,2 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.2')
testnr = 13
m1 = (/ 4 /)
z1 = dcmplx((/ 0 /))
y = 1.6
k = 1
res = GPL(m1,z1,y,k)
ref = dcmplx(0.0173042341866201179)
call check(res,ref)
call test_one_GPL((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3')
end subroutine do_GPL_tests
subroutine check_one_GPL(m,z,y,k,ref)
! checks one GPL(m,z,y,k) against reference value ref
integer :: m(:), k
complex(kind=prec) :: z(:), y, res, ref
res = GPL(m,z,y,k)
call check(res,ref)
end subroutine check_one_GPL
END PROGRAM TEST
\ No newline at end of file
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