Commit ce2bbf7a authored by Luca's avatar Luca

zero arguments and zeta function

parent 0e52465a
...@@ -12,6 +12,15 @@ CONTAINS ...@@ -12,6 +12,15 @@ CONTAINS
res = merge(1,n*factorial(n-1),n==1) res = merge(1,n*factorial(n-1),n==1)
END FUNCTION factorial END FUNCTION factorial
FUNCTION zeta(n)
real(kind=prec) :: values(9), zeta
integer :: n
values = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, &
1.03692775514337, 1.0173430619844488, 1.008349277381923, &
1.0040773561979441, 1.0020083928260821, 1.000994575127818/)
zeta = values(n-1)
END FUNCTION zeta
FUNCTION GPL_has_convergent_series(m,z,y,k) FUNCTION GPL_has_convergent_series(m,z,y,k)
! tests if GPL has a convergent series representation ! tests if GPL has a convergent series representation
integer :: m(:), k integer :: m(:), k
...@@ -34,6 +43,7 @@ CONTAINS ...@@ -34,6 +43,7 @@ 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
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
...@@ -48,7 +58,7 @@ CONTAINS ...@@ -48,7 +58,7 @@ CONTAINS
! are all z_i = 0 ? ! are all z_i = 0 ?
if(k == 1 .and. z(1) == 0) then if(k == 1 .and. z(1) == 0) then
! for that we assume that only one argument was passed, the rest through m1 ! for that we assume that only one argument was passed, the rest through m1
GPL = GPL_zero_zi(m(1)-1,y) GPL = GPL_zero_zi(m(1),y)
return return
end if end if
......
...@@ -6,21 +6,22 @@ PROGRAM TEST ...@@ -6,21 +6,22 @@ PROGRAM TEST
use mpl_module use mpl_module
use gpl_module use gpl_module
implicit none implicit none
complex(kind=prec) :: res complex(kind=prec) :: res
real, parameter :: tol = 1.0e-14 real, parameter :: tol = 1.0e-14
logical :: tests_successful = .true. logical :: tests_successful = .true.
call do_MPL_tests() ! call do_MPL_tests()
call do_GPL_tests() call do_GPL_tests()
if(tests_successful) then if(tests_successful) then
print*, 'All tests passed. ' print*, 'All tests passed. '
else else
print*, 'Some tests failed. ' print*, 'Some tests failed. '
stop 1 stop 1
end if end if
CONTAINS CONTAINS
subroutine check(res, ref) subroutine check(res, ref)
...@@ -40,22 +41,19 @@ CONTAINS ...@@ -40,22 +41,19 @@ CONTAINS
integer :: m(:) integer :: m(:)
complex(kind=prec) :: x(:), ref complex(kind=prec) :: x(:), ref
character(len=*) :: test_id character(len=*) :: test_id
print*, ' ', 'testing MPL ', test_id, ' ...' print*, ' ', 'testing MPL ', test_id, ' ...'
res = MPL(m,x) res = MPL(m,x)
call check(res,ref) call check(res,ref)
end subroutine test_one_MPL end subroutine test_one_MPL
subroutine do_MPL_tests() subroutine do_MPL_tests()
integer :: m2(2), m3(3) complex(kind=prec) :: ref
integer :: bla
complex(kind=prec) :: x2(2), x3(3)
complex(kind=prec) :: res, ref
print*, 'doing MPL tests...' print*, 'doing MPL tests...'
ref = dcmplx(0.022696600480693277651633) ref = dcmplx(0.022696600480693277651633)
call test_one_MPL((/ 1,1 /),cmplx((/ 0.3156498673740053, 0.3431255827785649 /)),ref, '1.1') call test_one_MPL((/ 1,1 /),cmplx((/ 0.3156498673740053, 0.3431255827785649 /)),ref, '1.1')
ref = dcmplx(0.00023134615630308335448329926098409) ref = dcmplx(0.00023134615630308335448329926098409)
call test_one_MPL((/ 1,1 /),cmplx((/ 0.03, 0.5012562893380046 /)),ref, '1.2') call test_one_MPL((/ 1,1 /),cmplx((/ 0.03, 0.5012562893380046 /)),ref, '1.2')
...@@ -75,19 +73,18 @@ CONTAINS ...@@ -75,19 +73,18 @@ CONTAINS
end subroutine test_one_GPL end subroutine test_one_GPL
subroutine do_GPL_tests() subroutine do_GPL_tests()
integer :: m2(2), m1(1), k complex(kind=prec) :: ref
complex(kind=prec) :: z2(2), z1(1), y, res, ref
print*, 'doing GPL tests...' print*, 'doing GPL tests...'
ref = dcmplx(0.0819393734128676) ref = dcmplx(0.0819393734128676)
call test_one_GPL((/ 1,1 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.1') call test_one_GPL((/ 1,1 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.1')
ref = dcmplx(0.01592795952537145) ref = dcmplx(0.01592795952537145)
call test_one_GPL((/ 3,2 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.2') call test_one_GPL((/ 3,2 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.2')
ref = dcmplx(0.0173042341866201179) 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')
end subroutine do_GPL_tests end subroutine do_GPL_tests
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!
Please register or to comment