gpl_module.f90 1.99 KB
Newer Older
Luca Naterop's avatar
Luca Naterop committed
1

Luca's avatar
Luca committed
2 3 4 5 6 7 8
MODULE gpl_module
  use mpl_module

  implicit none

CONTAINS 

Luca's avatar
Luca committed
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

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's avatar
Luca committed
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's avatar
Luca Naterop committed
41 42
  FUNCTION GPL(m,z,y,k)
    ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y)
Luca's avatar
Luca committed
43
    ! assumes zero arguments expressed through the m's
Luca Naterop's avatar
Luca Naterop committed
44 45 46 47

    integer :: m(:), k, i
    complex(kind=prec) :: z(:), x(k), y, GPL

48
    ! are all z_i = 0 ? 
Luca's avatar
Luca committed
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
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's avatar
Luca Naterop committed
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's avatar
Luca committed
64 65 66 67

  END FUNCTION GPL

END MODULE gpl_module
Luca Naterop's avatar
Luca Naterop committed
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's avatar
Luca Naterop committed
76
!   complex(kind=prec) :: res, ref
Luca Naterop's avatar
Luca Naterop committed
77 78
  
!   res = GPL(m,z,y,2)
Luca Naterop's avatar
Luca Naterop committed
79
!   ref = dcmplx(0.0819393734128676)
Luca Naterop's avatar
Luca Naterop committed
80
!   print*, 'res=',res
Luca Naterop's avatar
Luca Naterop committed
81
!   print*, 'ref=',ref
Luca Naterop's avatar
Luca Naterop committed
82 83 84

! END PROGRAM test