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

Luca's avatar
Luca committed
2 3
MODULE gpl_module
  use mpl_module
4
  use utils
Luca's avatar
Luca committed
5 6 7 8
  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

Luca's avatar
Luca committed
15 16 17 18 19 20 21 22 23
  FUNCTION zeta(n) 
    real(kind=prec) :: values(9), zeta
    integer :: n
    values = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, &
               1.03692775514337, 1.0173430619844488, 1.008349277381923, & 
               1.0040773561979441, 1.0020083928260821, 1.000994575127818/)
    zeta = values(n-1)
  END FUNCTION zeta

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
  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
41 42 43 44
  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
45
    print*, 'computed value using zi = 0'
Luca's avatar
Luca committed
46 47 48
    GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l
  END FUNCTION GPL_zero_zi

49
  RECURSIVE FUNCTION G_flat(z_flat,y) result(res)
50
    ! Calls G function with flat arguments, that is, zeroes not passed through the m's. 
51
    complex(kind=prec) :: z_flat(:), y, res
52 53
    complex(kind=prec), allocatable :: z(:), s(:,:)
    integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i
54
    integer, allocatable :: m(:)
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79

    print*, 'G_flat called with args', abs(z_flat)

    ! remove trailing zeroes
    k = size(z_flat)
    kminusj = find_amount_trailing_zeros(z_flat)
    j = k - kminusj
    if(all(abs(z_flat) < zero)) then 
      print*, 'all are zero', abs(z_flat)
      res = GPL_zero_zi(k,y)
      return
    else if(kminusj > 0) then
      print*, 'we have',kminusj,'trailing zeroes'
      allocate(s(j,j))
      s = shuffle_with_zero(z_flat(1:j-1))
      res = log(y)*G_flat(z_flat(1:size(z_flat)-1),y)
      do i = 1,size(s,1)
        res = res - G_flat([s(i,:),z_flat(j),zero_array(kminusj-1)], y)
      end do
      res = res / kminusj
      deallocate(s)
      return
    end if

    ! transform to condensed notation
80 81 82 83 84 85 86 87 88 89
    m_prime = get_condensed_m(z_flat)
    if(find_first_zero(m_prime) == -1) then
      condensed_size = size(m_prime)
    else
      condensed_size = find_first_zero(m_prime)-1 
    end if
    allocate(m(condensed_size))
    allocate(z(condensed_size))
    m = m_prime(1:condensed_size)
    z = get_condensed_z(m,z_flat)
90 91

    res = G_condensed(m,z,y,size(m))
92 93
    deallocate(m)
    deallocate(z)
94
  END  FUNCTION G_flat
95

96
  RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res)
Luca Naterop's avatar
Luca Naterop committed
97
    ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y)
Luca's avatar
Luca committed
98
    ! assumes zero arguments expressed through the m's
Luca's avatar
Luca committed
99
    
100
    integer :: m(:), k, i
101 102
    complex(kind=prec) :: z(:), x(k), y, res, c(sum(m)+1,sum(m)+1), z_flat(sum(m)), a(sum(m)-1)

103 104 105 106
    print*, 'called G_condensed with args'
    print*, 'm = ', m
    print*, 'z = ', abs(z)

107
    ! are all z_i = 0 ? 
108
    if(k == 1 .and. abs(z(1)) < zero) then
109 110
      ! assumes that the zeros at the beginning are passed through m1
      res = GPL_zero_zi(m(1),y)
Luca's avatar
Luca committed
111 112
      return
    end if
113

114 115 116
    ! has trailing zeroes?
    if(abs(z(k)) < zero ) then
      ! we remove them in flat form
117
      z_flat = get_flattened_z(m,z)
118
      res = G_flat(z_flat,y)
119 120
    end  if

Luca's avatar
Luca committed
121
    ! need make convergent?
122
    if(.not. GPL_has_convergent_series(m,z,y,k)) then
Luca's avatar
Luca committed
123 124 125 126 127 128
      print*, 'need to make convergent'
      if(k == 1 .and. m(1) == 1) then
        print*, 'now we use the easy case. ha. ha.'
      else
        print*, '  ', 'does not have convergent series representation'
      end if
129 130
      res = 0
      return
131 132
    end if

Luca Naterop's avatar
Luca Naterop committed
133 134 135
    do i = 1,k
      x(i) = merge(y/z(1), z(i-1)/z(i),i == 1)
    end do
136 137
    print*, 'computed using MPL'
    res = (-1)**k * MPL(m,x)
138
  END FUNCTION G_condensed
Luca's avatar
Luca committed
139 140

END MODULE gpl_module
Luca Naterop's avatar
Luca Naterop committed
141 142