mpl_module.f90 2.17 KB
 Luca committed Apr 04, 2019 1 2 `````` MODULE mpl_module `````` Luca committed Apr 03, 2019 3 `````` implicit none `````` Luca committed Apr 03, 2019 4 5 `````` integer, parameter :: prec = selected_real_kind(15,32) `````` Luca committed Apr 03, 2019 6 `````` integer :: GPLInfinity = 30 ! the default n if it is not passed `````` Luca committed Apr 03, 2019 7 8 9 10 11 12 `````` CONTAINS FUNCTION dilog(x,n) ! Computes the dilog Li_2(x) using the series representation up to order n integer :: n `````` Luca committed Apr 03, 2019 13 `````` complex(kind=prec) :: x, dilog `````` Luca committed Apr 03, 2019 14 15 16 17 18 19 20 `````` integer :: i integer :: j(n) j = (/(i, i=1,n,1)/) dilog = sum(x**j / j**2) END FUNCTION dilog `````` Luca committed Apr 03, 2019 21 `````` FUNCTION polylog(m,x,n) `````` Luca committed Apr 03, 2019 22 `````` ! Computes the classical polylogarithm Li_m(x) using series representation up to order n `````` Luca committed Apr 03, 2019 23 `````` integer :: m `````` Luca committed Apr 03, 2019 24 `````` complex(kind=prec) :: x, polylog `````` Luca committed Apr 03, 2019 25 26 `````` integer :: i,n integer, allocatable :: j(:) `````` Luca committed Apr 03, 2019 27 28 29 30 31 `````` j = (/(i, i=1,n,1)/) polylog = sum(x**j / j**m) END FUNCTION polylog `````` Luca committed Apr 09, 2019 32 `````` FUNCTION MPL_converges(m,x) `````` Luca committed Apr 03, 2019 33 34 `````` ! checks if the MPL converges complex(kind=prec) :: x(:) `````` Luca committed Apr 09, 2019 35 `````` integer :: m(:) `````` Luca committed Apr 04, 2019 36 `````` logical :: MPL_converges `````` Luca committed Apr 09, 2019 37 `````` MPL_converges = .false. `````` Luca committed Apr 03, 2019 38 `````` if(abs(product(x)) < 1) then `````` Luca committed Apr 09, 2019 39 40 41 `````` if(m(1) /= 1 .or. x(1) /= 1) then MPL_converges = .true. end if `````` Luca committed Apr 03, 2019 42 `````` end if `````` Luca committed Apr 04, 2019 43 `````` END FUNCTION MPL_converges `````` Luca committed Apr 03, 2019 44 `````` `````` Luca committed Apr 04, 2019 45 `````` recursive FUNCTION MPL(m, x, n_passed) result(res) `````` Luca committed Apr 03, 2019 46 47 `````` ! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n integer :: m(:) `````` Luca committed Apr 03, 2019 48 49 `````` complex(kind=prec) :: x(:) complex(kind=prec) :: res `````` Luca committed Apr 03, 2019 50 51 `````` integer, optional :: n_passed integer :: i, n `````` Luca committed Apr 03, 2019 52 `````` `````` Luca committed Apr 03, 2019 53 `````` n = merge(n_passed,GPLInfinity,present(n_passed)) `````` Luca committed Apr 04, 2019 54 `````` `````` Luca committed Apr 03, 2019 55 56 57 58 59 60 61 62 63 64 65 `````` if(size(m) /= size(x)) then print*, 'Error: m and x must have the same length' end if if(size(m) == 1) then ! base case res = polylog(m(1),x(1), n) else ! recursion step res = 0 do i = 1, n `````` Luca committed Apr 04, 2019 66 `````` res = res + x(1)**i / i**m(1) * MPL(m(2:), x(2:), i - 1) `````` Luca committed Apr 03, 2019 67 68 69 70 `````` end do ! a nicer way to do it would be but problem is i ! i = (/(j, j=1,n, 1)/) `````` Luca Naterop committed Apr 05, 2019 71 72 `````` ! res = sum( x(1)**i / i**m(1) * MPL(m(2:), x(2:), i(1) - 1) ) ! we could get around this problem by rewriting MPL to operate on each i and returning an array `````` Luca committed Apr 03, 2019 73 74 `````` end if `````` Luca committed Apr 04, 2019 75 `````` END FUNCTION MPL `````` Luca committed Apr 03, 2019 76 `````` `````` Luca committed Apr 04, 2019 77 ``````END MODULE mpl_module `````` Luca committed Apr 03, 2019 78 `````` `````` Luca committed Apr 03, 2019 79 ``````! PROGRAM test `````` Luca committed Apr 04, 2019 80 ``````! use mpl_module `````` Luca committed Apr 03, 2019 81 ``````! logical :: result `````` Luca committed Apr 09, 2019 82 ``````! result = MPL_converges( dcmplx((/1.0d0,.7d0,.3d0/)), (/ 1,2,1 /) ) `````` Luca committed Apr 03, 2019 83 84 ``````! print*, result ! end PROGRAM test``````