mpl_module.f90 1.25 KB
 Luca committed Apr 04, 2019 1 2 `````` MODULE mpl_module `````` Luca committed May 08, 2019 3 4 `````` use globals use utils `````` Luca Naterop committed Jul 08, 2019 5 `````` use maths_functions `````` Luca committed Apr 03, 2019 6 `````` implicit none `````` Luca committed Apr 03, 2019 7 `````` `````` Luca committed Apr 03, 2019 8 9 ``````CONTAINS `````` Luca committed Apr 09, 2019 10 `````` FUNCTION MPL_converges(m,x) `````` Luca committed Apr 03, 2019 11 12 `````` ! checks if the MPL converges complex(kind=prec) :: x(:) `````` Luca committed Apr 09, 2019 13 `````` integer :: m(:) `````` Luca committed Apr 04, 2019 14 `````` logical :: MPL_converges `````` Luca committed Apr 09, 2019 15 `````` MPL_converges = .false. `````` Luca committed Apr 03, 2019 16 `````` if(abs(product(x)) < 1) then `````` Luca Naterop committed Jul 08, 2019 17 `````` if(m(1) /= 1 .or. abs(x(1) - 1) < zero) then `````` Luca committed Apr 09, 2019 18 19 `````` MPL_converges = .true. end if `````` Luca committed Apr 03, 2019 20 `````` end if `````` Luca committed Apr 04, 2019 21 `````` END FUNCTION MPL_converges `````` Luca committed Apr 03, 2019 22 `````` `````` Luca committed Apr 04, 2019 23 `````` recursive FUNCTION MPL(m, x, n_passed) result(res) `````` Luca committed Apr 03, 2019 24 25 `````` ! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n integer :: m(:) `````` Luca committed Apr 03, 2019 26 27 `````` complex(kind=prec) :: x(:) complex(kind=prec) :: res `````` Luca committed Apr 03, 2019 28 29 `````` integer, optional :: n_passed integer :: i, n `````` Luca committed Apr 03, 2019 30 `````` `````` Luca Naterop committed Jul 10, 2019 31 `````` n = merge(n_passed,MPLInfinity,present(n_passed)) `````` Luca committed Apr 04, 2019 32 `````` `````` Luca committed Apr 03, 2019 33 34 35 `````` if(size(m) /= size(x)) then print*, 'Error: m and x must have the same length' end if `````` Luca Naterop committed Jul 09, 2019 36 37 38 39 `````` res = 0 `````` Luca committed Apr 03, 2019 40 `````` if(size(m) == 1) then `````` Luca Naterop committed Jul 09, 2019 41 42 43 44 45 46 `````` ! 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 `````` Luca committed Apr 03, 2019 47 48 49 `````` else ! recursion step do i = 1, n `````` Luca committed Apr 04, 2019 50 `````` res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1) `````` Luca committed Apr 03, 2019 51 52 53 `````` end do end if `````` Luca committed Apr 04, 2019 54 `````` END FUNCTION MPL `````` Luca committed Apr 03, 2019 55 `````` `````` Luca committed Apr 04, 2019 56 ``END MODULE mpl_module``