Commit aea74d41 authored by Luca Naterop's avatar Luca Naterop

polylogs cases where they are real

parent 84a917a7
......@@ -5,4 +5,5 @@
*.out
*.o
*.mod
*.smod
\ No newline at end of file
*.smod
.DS_Store
\ No newline at end of file
......@@ -9,7 +9,7 @@ PROGRAM eval
call parse_cmd_args()
res = GPL([0,1,3,2]) ! should be 0.09593041677639341 - 0.8829351795197851*I
res = GPL([2,1,3,4]) ! here's an example where we need to shuffle away from last place
print*, res
! res = pending_integral(cmplx([1,0]) + epsilon ,1, cmplx([3,2]) + epsilon)
......
......@@ -13,7 +13,7 @@ CONTAINS
complex(kind=prec) :: x, res
integer :: i,n
integer, allocatable :: j(:)
n = 300
n = 1000
j = (/(i, i=1,n,1)/)
res = sum(x**j / j**m)
END FUNCTION naive_polylog
......@@ -109,21 +109,19 @@ CONTAINS
Li2 = Real(LI2_OLD,prec)
END FUNCTION Li2
FUNCTION dilog_in_unit_circle(x) result(res)
! evaluates for any argument x in unit circle
complex(kind=prec) :: x, res
res = naive_polylog(2,x)
END FUNCTION dilog_in_unit_circle
FUNCTION dilog(x) result(res)
RECURSIVE FUNCTION dilog(x) result(res)
! evaluates dilog for any argument
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(x) <= 1.0) then
res = naive_polylog(2,x)
if(abs(aimag(x)) < zero ) then
res = Li2(real(x))
else
res = naive_polylog(2,x)
endif
else
res = -dilog_in_unit_circle(1/x) - pi**2/6 - log(add_ieps(-x))**2 / 2
res = -dilog(1/x) - pi**2/6 - log(add_ieps(-x))**2 / 2
end if
END FUNCTION dilog
......@@ -235,9 +233,12 @@ CONTAINS
! evaluates trilog for any argument
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(x) <= 1.0) then
res = naive_polylog(3,x)
if(abs(aimag(x)) < zero ) then
res = Li3(real(x))
else
res = naive_polylog(3,x)
endif
else
res = naive_polylog(3,sub_ieps(x)**(-1)) - (log(-sub_ieps(x)))**3/6 - pi**2/6 * log(-sub_ieps(x))
end if
......@@ -247,6 +248,9 @@ CONTAINS
! computes the polylog for now naively (except for dilog half-naively)
integer :: m
complex(kind=prec) :: x,res
print*, 'called polylog with m = ', m
if(verb >= 70) print*, 'called polylog(',m,',',x,')'
if(m == 2) then
res = dilog(x)
......
......@@ -11,7 +11,7 @@ PROGRAM TEST
implicit none
complex(kind=prec) :: res
real, parameter :: tol = 1.0e-14
real, parameter :: tol = 1.0e-10
logical :: tests_successful = .true.
integer :: i
......@@ -39,7 +39,7 @@ CONTAINS
if(delta < tol) then
print*, ' ',' passed with delta = ', delta
else
print*, ' ',' failed with delta = ', delta
print*, ' ',' FAILED with delta = ', delta
tests_successful = .false.
end if
end subroutine check
......@@ -78,12 +78,12 @@ CONTAINS
call check(res,ref)
end subroutine test_one_condensed
subroutine test_one_flat(z,y,ref,test_id)
complex(kind=prec) :: z(:), y, res, ref
subroutine test_one_flat(z,ref,test_id)
complex(kind=prec) :: z(:), res, ref
character(len=*) :: test_id
print*, ' ', 'testing GPL ', test_id, ' ...'
res = G_flat(z,y)
res = GPL(z)
call check(res,ref)
end subroutine test_one_flat
......@@ -101,6 +101,8 @@ CONTAINS
ref = dcmplx(0.0020332632172573974)
call test_one_condensed((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3')
ref = cmplx((0.09593041677639341, -0.8829351795197851))
call test_one_flat(cmplx([0,1,3,2]),ref,'2.4')
end subroutine do_GPL_tests
subroutine do_shuffle_tests()
......
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