gpl_module.f90 1.99 KB
 Luca Naterop committed Apr 04, 2019 1 `````` `````` Luca committed Apr 04, 2019 2 3 4 5 6 7 8 ``````MODULE gpl_module use mpl_module implicit none CONTAINS `````` Luca committed Apr 09, 2019 9 10 11 12 13 14 `````` RECURSIVE FUNCTION factorial(n) result(res) integer, intent(in) :: n integer :: res res = merge(1,n*factorial(n-1),n==1) END FUNCTION factorial `````` Luca committed Apr 09, 2019 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 `````` FUNCTION GPL_has_convergent_series(m,z,y,k) ! tests if GPL has a convergent series representation integer :: m(:), k complex(kind=prec) :: 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 = (y/z(1) /= 1) else GPL_has_convergent_series = .true. end if end if END FUNCTION GPL_has_convergent_series `````` Luca committed Apr 09, 2019 33 34 35 36 37 38 39 40 `````` FUNCTION GPL_zero_zi(l,y) ! used to compute the value of GPL when all zi are zero integer :: l complex(kind=prec) :: y, GPL_zero_zi GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l END FUNCTION GPL_zero_zi `````` Luca Naterop committed Apr 04, 2019 41 42 `````` FUNCTION GPL(m,z,y,k) ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) `````` Luca committed Apr 09, 2019 43 `````` ! assumes zero arguments expressed through the m's `````` Luca Naterop committed Apr 04, 2019 44 45 46 47 `````` integer :: m(:), k, i complex(kind=prec) :: z(:), x(k), y, GPL `````` Luca committed Apr 09, 2019 48 `````` ! are all z_i = 0 ? `````` Luca committed Apr 09, 2019 49 50 51 52 53 `````` if(k == 1 .and. z(1) == 0) then ! for that we assume that only one argument was passed, the rest through m1 GPL = GPL_zero_zi(m(1)-1,y) return end if `````` Luca committed Apr 09, 2019 54 55 56 57 58 59 `````` ! do they have convergent series rep? if(.not. GPL_has_convergent_series(m,z,y,k)) then print*, ' ', 'does not have convergent series representation' end if `````` Luca Naterop committed Apr 04, 2019 60 61 62 63 `````` do i = 1,k x(i) = merge(y/z(1), z(i-1)/z(i),i == 1) end do GPL = (-1)**k * MPL(m,x) `````` Luca committed Apr 04, 2019 64 65 66 67 `````` END FUNCTION GPL END MODULE gpl_module `````` Luca Naterop committed Apr 04, 2019 68 69 70 71 72 73 74 75 `````` ! PROGRAM test ! ! used to test this module ! use gpl_module ! integer :: m(2) = (/ 1,1 /) ! complex(kind=prec) :: z(2) = dcmplx((/ 1.3d0, 1.1d0 /)) ! complex(kind=prec) :: y = 0.4 `````` Luca Naterop committed Apr 05, 2019 76 ``````! complex(kind=prec) :: res, ref `````` Luca Naterop committed Apr 04, 2019 77 78 `````` ! res = GPL(m,z,y,2) `````` Luca Naterop committed Apr 05, 2019 79 ``````! ref = dcmplx(0.0819393734128676) `````` Luca Naterop committed Apr 04, 2019 80 ``````! print*, 'res=',res `````` Luca Naterop committed Apr 05, 2019 81 ``````! print*, 'ref=',ref `````` Luca Naterop committed Apr 04, 2019 82 83 84 `````` ! END PROGRAM test ``````