Commit 9ff6a342 authored by ulrich_y's avatar ulrich_y

New MPL routine

This is the algorithm by GiNaC. In theory one
could extend this to add a caching mechanism such as

    complex(kind=prec) :: cache(size(x),MPLMaxQ)

    do q=1,j
      cache(:,q) = x**q/q**m
    enddo
    do q=1,MPLMaxQ
      res = t(1)

      ! Fortran uses Column-major order, hence cache(:,q) is
      ! faster than cache(q,:).
      cache(:,q+j-1) = x**(q+j-1)/(q+j-1)**m
      t(j) = t(j) + cache(j,q)
      do k=1,j-1
        t(j-k) = t(j-k) + t(j-k+1) * cache(j-k,k+q)
      enddo

      if (mod(q,2) .eq. 1) then
        if (abs(t(1)-res).lt.MPLdel) exit
      endif
    enddo

In practice this doesn't really help because any
time saved with the cache is paid back through the
allocation and clearing of cache(:,:). Both
variations work similarly well now. If at some
point we might need MPLs with many more arguments
(size(x)), this might change.
parent ee4b3891
......@@ -8,7 +8,7 @@ MODULE globals
real, parameter :: pi = 3.14159265358979323846
! The following parameters control the accuracy of the evaluation
integer, protected :: MPLInfinity = 30 ! the default outermost expansion order for MPLs, formerly GPLInfinity
real(kind=prec), protected :: MPLdel = 1e-15 ! if the MPL sum changes less then del it is truncated.
integer, protected :: PolylogInfinity = 1000 ! expansion order for Polylogs
real(kind=prec), protected :: HoelderCircle = 1.1 ! when to apply Hoelder convolution?
......@@ -36,10 +36,10 @@ CONTAINS
END SUBROUTINE parse_cmd_args
#endif
SUBROUTINE SET_OPTIONS(mplinf, liinf, hcircle)
real(kind=prec), optional :: hcircle
integer, optional :: mplinf, liinf
if (present(mplinf)) MPLInfinity = mplinf
SUBROUTINE SET_OPTIONS(mpldel, liinf, hcircle)
real(kind=prec), optional :: hcircle, mpldel
integer, optional :: liinf
if (present(mpldel)) MPLdel = mpldel
if (present(liinf)) PolyLogInfinity = liinf
if (present(hcircle)) HoelderCircle = hcircle
END SUBROUTINE
......
......@@ -20,37 +20,36 @@ CONTAINS
end if
END FUNCTION MPL_converges
recursive FUNCTION MPL(m, x, n_passed) result(res)
! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
FUNCTION MPL(m, x) result(res)
integer :: m(:)
complex(kind=prec) :: x(:)
complex(kind=prec) :: res
integer, optional :: n_passed
integer :: i, n
complex(kind=prec) :: t(size(x))
integer :: q, j, k
j = size(x)
n = merge(n_passed,MPLInfinity,present(n_passed))
if(size(m) /= size(x)) then
print*, 'Error: m and x must have the same length'
end if
res=0
q=0
t=0
do while(.true.)
res = t(1)
q = q+1
t(j) = t(j) + x(j)**q/q**m(j)
do k=1,j-1
t(j-k) = t(j-k) + t(j-k+1) * x(j-k)**(k+q)/(k+q)**m(j-k)
enddo
if (mod(q,2) .eq. 1) then
if (abs(t(1)-res).lt.MPLdel) exit
endif
enddo
res = t(1)
res = 0
if(size(m) == 1) then
! base case, we get a polylog
do i = 1, n
if(i**m(1) < 0) return ! roll over
if(abs(x(1)**i) < 1.e-250) return
res = res + x(1)**i / i**m(1)
end do
else
! recursion step
do i = 1, n
res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1)
end do
end if
END FUNCTION MPL
END MODULE mpl_module
......@@ -82,6 +82,11 @@ CONTAINS
ref = cmplx(0.000023446106415452030937059124671151)
call test_one_MPL((/ 2,1,2 /),cmplx((/ 0.03, 0.5012562893380046, 55.3832 /)),ref, '1.3')
ref = cmplx(-0.06565799418838372)
call test_one_MPL((/1, 1/), cmplx((/-0.25,-2. /)), ref, '1.4')
ref = cmplx(-0.03199896396564833)
call test_one_MPL((/2, 1/), cmplx((/-0.25,-2. /)), ref, '1.4')
end subroutine do_MPL_tests
subroutine test_one_condensed(m,z,y,k,ref,test_id)
......
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