Commit d303027f authored by ulrich_y's avatar ulrich_y

Use successive difference for polylogs instead of fixed number of terms

parent 6bf49f69
...@@ -78,9 +78,17 @@ The compilation process creates the following results ...@@ -78,9 +78,17 @@ The compilation process creates the following results
difference between two successive terms at which the series difference between two successive terms at which the series
expansion (6) is truncated. expansion (6) is truncated.
* `integer LiInf = 1000` * `real(kind=prec) :: Lidel = 1e-15`
difference $`|\delta_i|`$ between two successive terms at which the
series expansion for
```math
{\rm Li}_n(x) = \sum_{i=1}^\infty \frac{x^i}{i^n}
= \sum_{i=1}^\infty \delta_i
```
is truncated.
number of terms in the expansion of classical polylogarithms.
* `real(kind=prec) :: hCircle = 1.1` * `real(kind=prec) :: hCircle = 1.1`
......
...@@ -456,6 +456,14 @@ CONTAINS ...@@ -456,6 +456,14 @@ CONTAINS
call test_one_flat( (/ (0._prec,0.), (1._prec,0.), (0.1_prec,-0.1_prec), (-1.3_prec,10.2_prec), (0.4_prec,0._prec) /) , ref, 'F.19') call test_one_flat( (/ (0._prec,0.), (1._prec,0.), (0.1_prec,-0.1_prec), (-1.3_prec,10.2_prec), (0.4_prec,0._prec) /) , ref, 'F.19')
ref = (0.6492441149301267712849185177612894261440132_prec,2.76692869057070821067245201199781913943_prec)
call iprint(' testing GPL G.1 ...',-1)
res = G([ toinum( (-1641._prec,4682._prec)/5000. ),&
toinum( (-1187._prec,4827._prec)/5000. ),&
toinum( 1069._prec/5000., +1_1 ), toinum( 787._prec / 5000., -1_1)], toinum(1._prec) )
call check(res,ref, ttol=5.e-9_prec)
end subroutine do_GPL_tests end subroutine do_GPL_tests
......
...@@ -16,7 +16,7 @@ MODULE globals ...@@ -16,7 +16,7 @@ MODULE globals
! The following parameters control the accuracy of the evaluation ! The following parameters control the accuracy of the evaluation
real(kind=prec), protected :: MPLdelta = zero ! if the MPL sum changes less then del it is truncated. real(kind=prec), protected :: MPLdelta = zero ! if the MPL sum changes less then del it is truncated.
integer, protected :: PolylogInfinity = 1000 ! expansion order for Polylogs real(kind=prec), protected :: Lidelta = zero ! like MPLdelta but for polylogs
real(kind=prec), protected :: HoelderCircle = 1.1_prec ! when to apply Hoelder convolution? real(kind=prec), protected :: HoelderCircle = 1.1_prec ! when to apply Hoelder convolution?
integer, parameter :: PolyLogCacheSize(2) = (/ 5, 100 /) integer, parameter :: PolyLogCacheSize(2) = (/ 5, 100 /)
! = (/ mmax, n /). At most n polylogs with weight mmax will be cached ! = (/ mmax, n /). At most n polylogs with weight mmax will be cached
...@@ -46,11 +46,10 @@ CONTAINS ...@@ -46,11 +46,10 @@ CONTAINS
END SUBROUTINE parse_cmd_args END SUBROUTINE parse_cmd_args
#endif #endif
SUBROUTINE SET_OPTIONS(mpldel, liinf, hcircle) SUBROUTINE SET_OPTIONS(mpldel, lidel, hcircle)
real(kind=prec), optional :: hcircle, mpldel real(kind=prec), optional :: hcircle, mpldel, lidel
integer, optional :: liinf
if (present(MPLdel)) MPLdelta = mpldel if (present(MPLdel)) MPLdelta = mpldel
if (present(liinf)) PolyLogInfinity = liinf if (present( Lidel)) LiDelta = lidel
if (present(hcircle)) HoelderCircle = hcircle if (present(hcircle)) HoelderCircle = hcircle
END SUBROUTINE END SUBROUTINE
......
...@@ -31,14 +31,19 @@ CONTAINS ...@@ -31,14 +31,19 @@ CONTAINS
FUNCTION naive_polylog(m,x) result(res) FUNCTION naive_polylog(m,x) result(res)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n ! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m integer :: m
complex(kind=prec) :: x, res complex(kind=prec) :: x, res, del
integer(kind=ikin) :: i integer(kind=ikin) :: i
res=0._prec res=0._prec
do i=1,PolylogInfinity
i = 1
del = 1._prec
do while (abs(del) > zero)
if(i**m.lt.0) return ! roll over if(i**m.lt.0) return ! roll over
if(abs(x**i).lt.1.e-250_prec) return if(abs(x**i).lt.1.e-250_prec) return
res = res+x**i/i**m del = x**i/i**m
enddo res = res+del
i = i+1
end do
END FUNCTION naive_polylog END FUNCTION naive_polylog
FUNCTION Li2(x) FUNCTION Li2(x)
......
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