Commit 014326e7 authored by ulrich_y's avatar ulrich_y

Merge branch 'mid-pending' into ieps

parents 468d2ea8 281ece3b
......@@ -12,9 +12,9 @@ PROGRAM eval
call parse_cmd_args()
res = GPL([0.1,-1.,1.])
res = GPL([-1.,0.,0.,0.,1.])
print*, res
! res = pending_integral(cmplx([1,3]),2,cmplx([2,4]))
! print*, res
......
......@@ -13,32 +13,6 @@ MODULE gpl_module
END INTERFACE GPL
CONTAINS
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)
! tests if GPL has a convergent series representation
integer :: m(:)
type(inum) :: 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 = .true. !(abs((y/z(1) - 1)) < zero)
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
......@@ -68,7 +42,7 @@ CONTAINS
if(.not. present(y)) print*, 'G(', abs(z_flat), ')'
END SUBROUTINE print_G
FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) result(res)
RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) result(res)
! here what is passed is not the full a vector, only a1, ..., ak without the trailing zeroes
integer :: m, i, j, n
type(inum) :: a(:), y2, s(m), p(:)
......@@ -133,7 +107,7 @@ CONTAINS
return
end if
a = g(1:size(p)-1)
a = g(1:size(g)-1)
y2 = g(size(g))
m = size(g) ! actually, size(g)-1+1
......@@ -187,11 +161,11 @@ CONTAINS
! case higher depth, s_r in middle, use my (67)
if(verb >= 30) print*, 's_r in the middle under PI'
res = -pending_integral(p,1,zeroes(0)) * G_flat([a(1:i-1),izero,a(i+1:size(a))],y2) &
+ pending_integral([p,a(i-1)],i-1,[a(1:i-2),a(i:size(a)),y2]) &
+ pending_integral([p,a(i-1)],1,zeroes(0)) * G_flat(a,y2) &
+ pending_integral([p,a(i)], i, [a(1:i-1), a(i+1:size(a)),y2]) &
- pending_integral([p,a(i)],1,zeroes(0)) * G_flat(a,y2)
res = +pending_integral(p,1,zeroes(0)) * G_flat([a(1:i-1),izero,a(i:size(a))],y2) & !64
- pending_integral([p,a(i-1)],i-1,[a(1:i-2),a(i:size(a)),y2]) & ! 67a
+ pending_integral([p,a(i-1)],1,zeroes(0)) * G_flat(a,y2) & ! 67c
+ pending_integral([p,a(i)], i, [a(1:i-1), a(i+1:size(a)),y2]) & !67b
- pending_integral([p,a(i)],1,zeroes(0)) * G_flat(a,y2) !67a
END FUNCTION pending_integral
RECURSIVE FUNCTION remove_sr_from_last_place_in_G(a,y2,m,sr) result(res)
......@@ -298,10 +272,9 @@ CONTAINS
if(verb >= 50) call print_G(z_flat,y)
! catch G(1,1)
if(size(z_flat) == 1) then
if( abs(z_flat(1) - y) <= zero ) then
print*, 'catch G(1,1)'
res = 0
return
end if
......@@ -440,7 +413,7 @@ CONTAINS
end if
! need make convergent?
if(.not. GPL_has_convergent_series(m,z,y)) then
if(.not. all(abs(y) <= abs(z))) then
z_flat = get_flattened_z(m,z)
res = G_flat(z_flat,y)
return
......
......@@ -8,16 +8,28 @@ MODULE maths_functions
module procedure polylogcmplx,polyloginum
end interface polylog
CONTAINS
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 naive_polylog(m,x) result(res)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m
complex(kind=prec) :: x, res
integer :: i,n
integer, allocatable :: j(:)
n = 1000
j = (/(i, i=1,n,1)/)
res = sum(x**j / j**m)
res=0.
do i=1,n
if(i**m.lt.0) return ! roll over
if(abs(x**i).lt.1.e-250) return
res = res+x**i/i**m
enddo
END FUNCTION naive_polylog
FUNCTION Li2(x)
......@@ -253,6 +265,48 @@ CONTAINS
res = polylog(m,x%c)
END FUNCTION polyloginum
FUNCTION BERNOULLI_POLYNOMIAL(n, x) result(res)
integer n
complex(kind=prec) :: x, res
select case(n)
case(1)
res = -1/2. + x
case(2)
res = 1/6. - x + x**2
case(3)
res = x/2. - 3*x**2/2. + x**3
case(4)
res = -1/30. + x**2 - 2*x**3 + x**4
case(5)
res = -x/6. + 5*x**3/3 - 5*x**4/2 + x**5
case(6)
res = 1/42. - x**2/2 + 5*x**4/2 - 3*x**5 + x**6
case(7)
res = x/6. - 7*x**3/6 + 7*x**5/2 - 7*x**6/2 + x**7
case(8)
res = -1/30. + 2*x**2/3 - 7*x**4/3 + 14*x**6/3 - 4*x**7 + x**8
case(9)
res = -3*x/10 + 2*x**3 - 21*x**5/5 + 6*x**7 - 9*x**8/2 + x**9
case(10)
res = 5/66. - 3*x**2/2 + 5*x**4 - 7*x**6 + 15*x**8/2 - 5*x**9 + x**10
case(11)
res = 5*x/6 - 11*x**3/2 + 11*x**5 - 11*x**7 + 55*x**9/6 - 11*x**10/2 + x**11
case(12)
res = -691/2730. + 5*x**2 - 33*x**4/2 + 22*x**6 - 33*x**8/2 + 11*x**10 - 6*x**11 + x**12
case(13)
res = -691*x/210 + 65*x**3/3 - 429*x**5/10 + 286*x**7/7 - 143*x**9/6 + 13*x**11 - 13*x**12/2 + x**13
case(14)
res = 7/6. - 691*x**2/30 + 455*x**4/6 - 1001*x**6/10 + 143*x**8/2 - 1001*x**10/30 + 91*x**12/6 - 7*x**13 + x**14
case(15)
res = 35*x/2 - 691*x**3/6 + 455*x**5/2 - 429*x**7/2 + 715*x**9/6 - 91*x**11/2 + 35*x**13/2 - 15*x**14/2 + x**15
case default
print*,"Bernoulli beyond 15 is not implemented"
stop
end select
END FUNCTION
FUNCTION polylogcmplx(m,x) result(res)
! computes the polylog
......@@ -260,12 +314,21 @@ CONTAINS
complex(kind=prec) :: x,res
if(verb >= 70) print*, 'called polylog(',m,',',x,')'
if(m == 2) then
if ((m.le.9).and.(abs(x-1.).lt.zero)) then
res = zeta(m)
else if ((m.le.9).and.(abs(x+1.).lt.zero)) then
res = -(1. - 2.**(1-m))*zeta(m)
else if(m == 2) then
res = dilog(x)
else if(m == 3) then
res = trilog(x)
else
res = naive_polylog(m,x)
if (abs(x).gt.1) then
res = (-1)**(m-1)*naive_polylog(m,1./x) &
- cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*clog(-x)/2/pi) / factorial(m)
else
res = naive_polylog(m,x)
endif
end if
END FUNCTION polylogcmplx
......
......@@ -2,6 +2,7 @@
MODULE mpl_module
use globals
use utils
use maths_functions
implicit none
CONTAINS
......
......@@ -130,8 +130,8 @@ CONTAINS
ref = (-8.403785974849544e-6,0.0)
call test_one_flat(cmplx([-1./xchen,0.,-1./xchen,1./(xchen*z),1.0]),ref,'3.5')
! ref = cmplx((0.4925755847450199,2.6389214054743295))
! call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
ref = cmplx((0.4925755847450199,2.6389214054743295))
call test_one_flat(cmplx([-1.,-1.,z,z,1.]),ref,'3.6')
! ref = cmplx((-0.0015317713178859967,-0.00045003911367000565))
! call test_one_flat(cmplx([0.,0.,(1-sqrt(1-z**2))/(xchen*z), 1./(xchen*z),1.]),ref,'3.7')
......@@ -177,7 +177,6 @@ CONTAINS
end subroutine do_GPL_tests
! subroutine do_shuffle_tests()
! complex(kind=prec) :: v(2) = cmplx((/1,2/))
! complex(kind=prec) :: w(2) = cmplx((/3,4/))
......
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