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

Luca's avatar
Luca committed
2
MODULE gpl_module
3
  use globals
ulrich_y's avatar
ulrich_y committed
4
  use ieps
ulrich_y's avatar
ulrich_y committed
5 6 7 8
  use utils
  use shuffle, only: shuffle_product, shuffle_with_zero
  use maths_functions, only: plog1, polylog, zeta
  use mpl_module, only: mpl
Luca's avatar
Luca committed
9 10
  implicit none

ulrich_y's avatar
ulrich_y committed
11
  INTERFACE G
ulrich_y's avatar
ulrich_y committed
12 13 14
    module procedure G_flat, G_flatR, G_flatC, &
                     G_condensed, G_condensedS, G_condensedR, G_condensedC, &
                     G_superflat, G_real, G_int
ulrich_y's avatar
ulrich_y committed
15
  END INTERFACE G
Luca's avatar
Luca committed
16 17
CONTAINS 

Luca's avatar
Luca committed
18 19 20
  FUNCTION GPL_zero_zi(l,y)
    ! used to compute the value of GPL when all zi are zero
    integer :: l
ulrich_y's avatar
ulrich_y committed
21 22
    type(inum) :: y
    complex(kind=prec) :: GPL_zero_zi
ulrich_y's avatar
ulrich_y committed
23 24 25 26
    if (abs(aimag(y)).lt.zero) then
      if (real(y).gt.0) then
        GPL_zero_zi = 1.0_prec/factorial(l) * log(real(y)) ** l
      else
ulrich_y's avatar
ulrich_y committed
27
        GPL_zero_zi = 1.0_prec/factorial(l) * (log(-real(y))+i_*(y%i0*pi)) ** l
ulrich_y's avatar
ulrich_y committed
28 29 30 31
      endif
    else
      GPL_zero_zi = 1.0_prec/factorial(l) * log(y%c) ** l
    endif
Luca's avatar
Luca committed
32 33
  END FUNCTION GPL_zero_zi

Luca's avatar
Luca committed
34 35
  FUNCTION is_convergent(z,y)
    ! returns true if G(z,y) convergent, otherwise false
Luca's avatar
Luca committed
36
    ! can be used in either notation (flat or condensed)
ulrich_y's avatar
ulrich_y committed
37
    type(inum) :: z(:), y
Luca's avatar
Luca committed
38 39 40 41 42 43 44 45 46 47
    logical :: is_convergent
    integer :: i

    is_convergent = .true.
    do i = 1,size(z)
      if(abs(z(i)) < zero) cycle  ! skip zero values
      if(abs(y) > abs(z(i))) is_convergent = .false.
    end do
  END FUNCTION is_convergent

Luca's avatar
Luca committed
48
  SUBROUTINE print_G(z_flat, y)
ulrich_y's avatar
ulrich_y committed
49 50
    type(inum)  :: z_flat(:)
    type(inum) , optional :: y
51
    if(present(y)) print*, 'G(', real(z_flat), real(y), ')'
52
    if(.not. present(y)) print*, 'G(', abs(z_flat), ')'
Luca's avatar
Luca committed
53 54
  END SUBROUTINE print_G

ulrich_y's avatar
ulrich_y committed
55
  RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p, srs) result(res)
56 57
    ! here what is passed is not the full a vector, only a1, ..., ak without the trailing zeroes
    integer :: m, i, j, n
ulrich_y's avatar
ulrich_y committed
58
    integer(1) :: srs
ulrich_y's avatar
ulrich_y committed
59 60
    type(inum) :: a(:), y2, s(m), p(:)
    complex(kind=prec) :: res
61
    type(inum) :: alpha(product((/(i,i=1,size(a)+size(s))/))/  &
62 63 64
        (product((/(i,i=1,size(a))/))*product((/(i,i=1,size(s))/))), & 
        size(a) + size(s))

ulrich_y's avatar
ulrich_y committed
65
    s = [zeroes(m-1),marker]
66
    alpha = shuffle_product(a,s)
ulrich_y's avatar
ulrich_y committed
67
#ifdef DEBUG
68 69 70
    if(verb >= 50) then
      print*, 'mapping to '
      call print_G(a,y2)
71
      print*, 'PI with p=',real(p),'i=',m,'g =',real([zeroes(m-1),y2])
72
    end if
ulrich_y's avatar
ulrich_y committed
73
#endif
ulrich_y's avatar
ulrich_y committed
74
    res = G_flat(a,y2)*pending_integral(p,m,[zeroes(m-1),y2], srs)
ulrich_y's avatar
ulrich_y committed
75
#ifdef DEBUG
76
    if(verb >= 50) print*, 'also mapping to'
ulrich_y's avatar
ulrich_y committed
77
#endif
78 79
    do j = 2,size(alpha, 1)
      ! find location of s_r
ulrich_y's avatar
ulrich_y committed
80
      n = find_marker(alpha(j,:))
ulrich_y's avatar
ulrich_y committed
81
#ifdef DEBUG
82 83
      if(verb >= 50) print*, 'PI with p=',real(p),'i=',n,'g =',&
          real([alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2])
ulrich_y's avatar
ulrich_y committed
84
#endif
ulrich_y's avatar
ulrich_y committed
85
      res = res - pending_integral(p, n, [alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2], srs)
86 87 88
    end do
  END FUNCTION remove_sr_from_last_place_in_PI

ulrich_y's avatar
ulrich_y committed
89
  RECURSIVE FUNCTION pending_integral(p,i,g,srs) result(res)
90
    ! evaluates a pending integral by reducing it to simpler ones and g functions
ulrich_y's avatar
ulrich_y committed
91 92 93
    complex(kind=prec) :: res
    type(inum) :: p(:), g(:)
    type(inum) :: y1, y2, b(size(p)-1), a(size(g)-1)
Luca's avatar
Luca committed
94
    integer :: i, m
ulrich_y's avatar
ulrich_y committed
95
    integer(1):: srs
96
    res = 0
97

ulrich_y's avatar
ulrich_y committed
98
#ifdef DEBUG
99
    if(verb >= 30) print*, 'evaluating PI with p=',real(p),'i=',real(i),'g =',real(g)
ulrich_y's avatar
ulrich_y committed
100
#endif
101 102 103 104 105 106

    y1 = p(1)
    b = p(2:size(p))

    ! if integration variable is not in G-function
    if(i == 0 .or. size(g) == 0) then
ulrich_y's avatar
ulrich_y committed
107
#ifdef DEBUG
108
      if(verb >= 30) print*, 'only integrals in front, we get G-function'
ulrich_y's avatar
ulrich_y committed
109
#endif
110 111 112
      res = G_flat(b,y1)
      return
    end if
Luca's avatar
Luca committed
113

114
    ! if integration variable at end -> we gat a G function 
115
    if(i == size(g)+1) then
ulrich_y's avatar
ulrich_y committed
116
#ifdef DEBUG
117
      if(verb >= 30) print*, 'is just a G-function'
ulrich_y's avatar
ulrich_y committed
118
#endif
119 120 121
      res = G_flat([p(2:size(p)),g], p(1))
      return
    end if
Luca's avatar
Luca committed
122
  
123

124
    ! if depth one and m = 1 use (23) 
125
    if(size(g) == 1) then
ulrich_y's avatar
ulrich_y committed
126
#ifdef DEBUG
127
      if(verb >= 30) print*, 'case depth one and m = 1'
ulrich_y's avatar
ulrich_y committed
128
#endif
ulrich_y's avatar
ulrich_y committed
129 130
      !res = pending_integral(p,2,[sub_ieps(g(1))]) - pending_integral(p,2,[cmplx(0.0)]) &
      !  + G_flat(p(2:size(p)), p(1)) * log(-sub_ieps(g(1)))
ulrich_y's avatar
ulrich_y committed
131 132
      res = pending_integral(p,2,[inum( g(1)%c,-srs )  ], srs) - pending_integral(p,2,[izero], srs) &
       + G_flat(p(2:size(p)), p(1)) * (log(g(1)%c) +  srs    * pi * i_)
133 134
      return
    end if
Luca's avatar
Luca committed
135
  
136
    a = g(1:size(g)-1)
Luca Naterop's avatar
Luca Naterop committed
137 138 139
    y2 = g(size(g)) 
    m = size(g)  ! actually, size(g)-1+1

140
    ! if depth one and m > 1 use (25)
141
    if(all( abs( g(1:size(g)-1) ) < zero)) then       
ulrich_y's avatar
ulrich_y committed
142
#ifdef DEBUG
143 144 145 146
      if(verb >= 30) print*, 'case depth one and m > 1'
      if(verb >= 50) then 
        print*, 'map to'
        print*, 'zeta(',m,')'
147 148 149
        print*, 'PI with p=',real(p),'i=',0,'g =',tocmplx(zeroes(0))
        print*, 'PI with p=',real([y2,izero]),'i=',m-1,'g =',tocmplx([zeroes(m-2),y2])
        print*, 'PI with p=',real(p),'i=',0,'g =',tocmplx(zeroes(0))
ulrich_y's avatar
ulrich_y committed
150
        print*, 'PI with p=',tocmplx([p,izero]),'i=',m-1,'g =',tocmplx(zeroes(0))
151
      end if
ulrich_y's avatar
ulrich_y committed
152
#endif
ulrich_y's avatar
ulrich_y committed
153 154 155
      res = -zeta(m)*pending_integral(p,0,zeroes(0),srs) &
        + pending_integral([y2,izero],m-1,[zeroes(m-2),y2],srs)*pending_integral(p,0,zeroes(0),srs) &
        - pending_integral([p,izero] ,m-1,[zeroes(m-2),y2],srs)
Luca's avatar
Luca committed
156 157 158
      return
    end if
  
159
    ! case of higher depth, s_r at beginning, use (21)
160
    if(i == 1) then
ulrich_y's avatar
ulrich_y committed
161
#ifdef DEBUG
162
      if(verb >= 30) print*, 'case higher depth, sr at beginning'
Luca's avatar
Luca committed
163
      
164
      if(verb >= 50) then
Luca Naterop's avatar
minor  
Luca Naterop committed
165
        print*, 'map to (using 21)'
166
        print*, 'PI with p=',real(p),'i=',0,'g =',tocmplx(zeroes(0) )
ulrich_y's avatar
ulrich_y committed
167
        call print_G([izero,a],y2)
168
        print*, 'PI with p=',real([p,y2]),'i=',0,'g =',tocmplx(zeroes(0) )
169
        call print_G(a,y2)
170
        print*, 'PI with p=',real([p,a(1)]),'i=',1,'g =',real([a(2:size(a)),y2])
171
        print*, 'PI with p=',real([p,a(1)]),'i=',0,'g =',tocmplx(zeroes(0))
172 173
        call print_G(a,y2)
      end if
ulrich_y's avatar
ulrich_y committed
174
#endif
175

ulrich_y's avatar
ulrich_y committed
176 177 178 179
      res = pending_integral(p,0,zeroes(0),srs) * G_flat([izero,a],y2) &
        + pending_integral([p,y2],0,zeroes(0),srs) * G_flat(a,y2) &
        + pending_integral([p,a(1)],1,[a(2:size(a)),y2],srs) &
        - pending_integral([p,a(1)],0,zeroes(0),srs) * G_flat(a,y2)
180 181
        return
    end if
Luca's avatar
Luca committed
182
    
183
    ! case higher depth, s_r at the end, use shuffle algebra to remove from last place
Luca's avatar
Luca committed
184
    if(i == size(g)) then
ulrich_y's avatar
ulrich_y committed
185
#ifdef DEBUG
Luca Naterop's avatar
Luca Naterop committed
186
      if(verb >= 30) print*, 's_r at the end under PI, need to shuffle'
ulrich_y's avatar
ulrich_y committed
187
#endif
188
      m = find_amount_trailing_zeros(a) + 1
ulrich_y's avatar
ulrich_y committed
189
      res = remove_sr_from_last_place_in_PI(a(1:size(a)-(m-1)), y2, m, p, srs)
Luca's avatar
Luca committed
190 191 192
      return
    end if
    
193
    ! case higher depth, s_r in middle, use (22)
ulrich_y's avatar
ulrich_y committed
194
#ifdef DEBUG
Luca Naterop's avatar
Luca Naterop committed
195
    if(verb >= 30) print*, 's_r in the middle under PI'
ulrich_y's avatar
ulrich_y committed
196
#endif
Luca Naterop's avatar
Luca Naterop committed
197

ulrich_y's avatar
ulrich_y committed
198 199 200 201 202
    res =  +pending_integral(p,1,zeroes(0),srs) * G_flat([a(1:i-1),izero,a(i:size(a))],y2) &
      - pending_integral([p,a(i-1)],i-1,[a(1:i-2),a(i:size(a)),y2],srs) &
      + pending_integral([p,a(i-1)],1,zeroes(0),srs) * G_flat(a,y2) &
      + pending_integral([p,a(i)], i, [a(1:i-1), a(i+1:size(a)),y2],srs)  &
      - pending_integral([p,a(i)],1,zeroes(0),srs) * G_flat(a,y2)
Luca's avatar
Luca committed
203 204
  END FUNCTION pending_integral

205
  RECURSIVE FUNCTION remove_sr_from_last_place_in_G(a,y2,m,sr) result(res)
ulrich_y's avatar
ulrich_y committed
206 207
    type(inum) :: a(:), sr, y2
    complex(kind=prec) :: res
Luca's avatar
Luca committed
208
    integer :: m,i,j
209
    type(inum) :: alpha(product((/(i,i=1,size(a)+m)/))/  &
Luca's avatar
Luca committed
210 211
      (product((/(i,i=1,size(a))/))*product((/(i,i=1,m)/))), & 
      size(a) + m)
Luca's avatar
Luca committed
212 213
    alpha = shuffle_product(a,[zeroes(m-1),sr])
    res = G_flat(a,y2)*G_flat([zeroes(m-1),sr],y2)
Luca's avatar
Luca committed
214 215 216
    do j = 2,size(alpha,1)
      res = res - G_flat(alpha(j,:),y2)
    end do
217
  END FUNCTION remove_sr_from_last_place_in_G
Luca's avatar
Luca committed
218

219
  RECURSIVE FUNCTION make_convergent(a,y2) result(res)
220
    ! goes from G-functions to pending integrals and simpler expressions
Luca's avatar
Luca committed
221

ulrich_y's avatar
ulrich_y committed
222 223
    type(inum) :: a(:), y2, sr
    complex(kind=prec) :: res
Luca's avatar
Luca committed
224
    integer :: i, mminus1
Luca's avatar
Luca committed
225

Luca's avatar
Luca committed
226
    res = 0
227 228
    i = min_index(abs(a))
    sr = a(i)
Luca's avatar
Luca committed
229

Luca Naterop's avatar
Luca Naterop committed
230
    if(i == size(a)) then
231
      ! sr at the end, thus shuffle
ulrich_y's avatar
ulrich_y committed
232
#ifdef DEBUG
Luca Naterop's avatar
Luca Naterop committed
233
      if(verb >= 30) print*, 'sr at the end'
ulrich_y's avatar
ulrich_y committed
234
#endif
Luca Naterop's avatar
Luca Naterop committed
235 236 237 238
      mminus1 = find_amount_trailing_zeros(a(1:size(a)-1))
      res = remove_sr_from_last_place_in_G(a(1:size(a)-mminus1-1),y2,mminus1+1,sr)
      return
    end if    
Luca's avatar
Luca committed
239

240
    if(i == 1) then
241
      !s_r at beginning, thus use (18)
Luca's avatar
Luca committed
242

ulrich_y's avatar
ulrich_y committed
243
#ifdef DEBUG
Luca's avatar
Luca committed
244 245 246
      if(verb >= 30) then 
        print*, '--------------------------------------------------'
        print*, 'sr at beginning, map to: '
ulrich_y's avatar
ulrich_y committed
247
        call print_G([izero, a(i+1:size(a))], y2)
Luca's avatar
Luca committed
248 249
        call print_G([y2], sr)
        call print_G(a(i+1:size(a)), y2)
250
        print*, 'PI with p=',real([sr, a(i+1)]),'i=',i,'g =', real([a(i+2:size(a)), y2])
Luca's avatar
Luca committed
251 252 253 254
        call print_G([a(i+1)], sr)
        call print_G(a(i+1:size(a)), y2)
        print*, '--------------------------------------------------'
      end if
ulrich_y's avatar
ulrich_y committed
255
#endif
ulrich_y's avatar
ulrich_y committed
256
      res = G_flat([izero, a(i+1:size(a))], y2) &
257
        + G_flat([y2], sr) * G_flat(a(i+1:size(a)), y2) &
ulrich_y's avatar
ulrich_y committed
258
        + pending_integral([sr, a(i+1)], i, [a(i+2:size(a)), y2],sr%i0) &
ulrich_y's avatar
ulrich_y committed
259
        - G_flat([a(i+1)],sr) * G_flat(a(i+1:size(a)), y2)
Luca's avatar
Luca committed
260 261
      return
    end if
Luca's avatar
Luca committed
262

263
    ! so s_r in middle, use (19)
ulrich_y's avatar
ulrich_y committed
264
#ifdef DEBUG
265 266 267
    if(verb >= 30) then 
      print*, '--------------------------------------------------'
      print*, 'sr in the middle, map to: '
ulrich_y's avatar
ulrich_y committed
268
      call print_G([a(1:i-1),izero,a(i+1:size(a))],y2)
269
      print*, 'PI with p=',real([sr,a(i-1)]),'i=', i-1,'g =', real([a(1:i-2),a(i+1:size(a)),y2])
270 271
      call print_G([a(i-1)],sr)
      call print_G([a(1:i-1),a(i+1:size(a))],y2)
272
      print*, 'and PI with p=',real([sr,a(i+1)]),'i=',i,'g =', real([a(1:i-1),a(i+2:size(a)),y2])
273 274 275 276
      call print_G([a(i+1)],sr)
      call print_G([a(1:i-1),a(i+1:size(a))],y2)
      print*, '--------------------------------------------------'
    end if
ulrich_y's avatar
ulrich_y committed
277
#endif
278

ulrich_y's avatar
ulrich_y committed
279
    res = G_flat([a(1:i-1),izero,a(i+1:size(a))],y2) &
ulrich_y's avatar
ulrich_y committed
280
      - pending_integral([sr,a(i-1)], i-1, [a(1:i-2),a(i+1:size(a)),y2],sr%i0)  &
Luca's avatar
Luca committed
281
      + G_flat([a(i-1)],sr) * G_flat([a(1:i-1),a(i+1:size(a))],y2)        &
ulrich_y's avatar
ulrich_y committed
282
      + pending_integral([sr,a(i+1)], i, [a(1:i-1),a(i+2:size(a)),y2],sr%i0)    &
Luca's avatar
Luca committed
283
      - G_flat([a(i+1)],sr) * G_flat([a(1:i-1),a(i+1:size(a))],y2)        
Luca's avatar
Luca committed
284
  END FUNCTION make_convergent
Luca's avatar
Luca committed
285

ulrich_y's avatar
ulrich_y committed
286
  RECURSIVE FUNCTION improve_convergence(z) result(res)
Luca Naterop's avatar
Luca Naterop committed
287
    ! improves the convergence by applying the Hoelder convolution to G(z1,...zk,1)
ulrich_y's avatar
ulrich_y committed
288 289
    type(inum) :: z(:),oneminusz(size(z))
    complex(kind=prec) :: res
Luca Naterop's avatar
Luca Naterop committed
290 291
    complex(kind=prec), parameter :: p = 2.0
    integer :: k, j
ulrich_y's avatar
ulrich_y committed
292
#ifdef DEBUG
Luca Naterop's avatar
Luca Naterop committed
293
    if(verb >= 30) print*, 'requires Hoelder convolution'
ulrich_y's avatar
ulrich_y committed
294
#endif
ulrich_y's avatar
ulrich_y committed
295 296 297 298
    ! In the Hoelder expression, all the (1-z) are -i0.. GiNaC does something
    ! different (and confusing, l. 1035). As we do, they usually would set
    ! i0 -> -z%i0. However, if Im[z] == 0 and Re[z] >= 1, they just set it to
    ! i0 -> +i0, be damned what it was before.
ulrich_y's avatar
ulrich_y committed
299
    do j=1,size(z)
ulrich_y's avatar
ulrich_y committed
300 301 302 303 304
      if ( (abs(aimag(z(j))) .lt. zero).and.( real(z(j)) .ge. 1) ) then
        oneminusz(j) = inum(1.-z(j)%c, +1)
      else
        oneminusz(j) = inum(1.-z(j)%c,-z(j)%i0)
      endif
ulrich_y's avatar
ulrich_y committed
305
    enddo
Luca Naterop's avatar
Luca Naterop committed
306 307
    k = size(z)

ulrich_y's avatar
ulrich_y committed
308 309
    res = G_flat(z,inum(1./p,di0)) ! first term of the sum
    res = res + (-1)**k * G_flat(oneminusz(k:1:-1), inum(1.-1/p,di0))
Luca Naterop's avatar
Luca Naterop committed
310
    do j = 1,k-1
ulrich_y's avatar
ulrich_y committed
311
      res = res + (-1)**j * G_flat(oneminusz(j:1:-1),inum(1.-1/p,di0)) * G_flat(z(j+1:k),inum(1./p,di0))
Luca Naterop's avatar
Luca Naterop committed
312 313 314
    end do
  END FUNCTION improve_convergence

315
  RECURSIVE FUNCTION G_flat(z_flat,y) result(res)
316
    ! Calls G function with flat arguments, that is, zeroes not passed through the m's. 
ulrich_y's avatar
ulrich_y committed
317
    type(inum) :: z_flat(:), y, znorm(size(z_flat))
ulrich_y's avatar
ulrich_y committed
318 319
    complex(kind=prec) :: res
    type(inum), allocatable :: s(:,:), z(:)
Luca's avatar
Luca committed
320
    integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i, m_1
321
    integer, allocatable :: m(:)
Luca's avatar
Luca committed
322
    logical :: is_depth_one
323

ulrich_y's avatar
ulrich_y committed
324
#ifdef DEBUG
325
    if(verb >= 50) call print_G(z_flat,y)
ulrich_y's avatar
ulrich_y committed
326
#endif
Luca's avatar
Luca committed
327

328

329
    if(size(z_flat) == 1) then
ulrich_y's avatar
ulrich_y committed
330
      if( abs(z_flat(1)%c - y%c) <= zero ) then
331 332 333 334 335
        res = 0
        return
      end if
    end if

336
    ! add small imaginary part if not there
Luca Naterop's avatar
Luca Naterop committed
337 338 339 340
    ! do i = 1,size(z_flat)
    !   if(abs(aimag(z_flat(i))) < 1e-25) z_flat(i) = add_ieps(z_flat(i))
    !   if(abs(aimag(y)) < 1e-25) y = add_ieps(y)
    ! end do
Luca's avatar
Luca committed
341

342
    ! is just a logarithm? 
Luca's avatar
Luca committed
343
    if(all(abs(z_flat) < zero)) then
ulrich_y's avatar
ulrich_y committed
344
#ifdef DEBUG
345
      if(verb >= 70) print*, 'all z are zero'
ulrich_y's avatar
ulrich_y committed
346
#endif
ulrich_y's avatar
ulrich_y committed
347
      res = gpl_zero_zi(size(z_flat),y)
Luca's avatar
Luca committed
348 349
      return
    end if
350
    if(size(z_flat) == 1) then
ulrich_y's avatar
ulrich_y committed
351
#ifdef DEBUG
352
      if(verb >= 70) print*, 'is just a logarithm'
ulrich_y's avatar
ulrich_y committed
353
#endif
ulrich_y's avatar
ulrich_y committed
354
      res = plog1(y,z_flat(1))  ! log((z_flat(1) - y)/z_flat(1))
355 356 357 358 359 360 361 362 363
      return
    end if

    ! is it a polylog? 
    m_prime = get_condensed_m(z_flat)
    m_1 = m_prime(1)
    is_depth_one = (count((m_prime>0)) == 1)
    if(is_depth_one) then
      ! case m >= 2, other already handled above
ulrich_y's avatar
ulrich_y committed
364
#ifdef DEBUG
Luca's avatar
Luca committed
365
      if(verb >= 70) print*, 'is just a polylog'
ulrich_y's avatar
ulrich_y committed
366
#endif
ulrich_y's avatar
ulrich_y committed
367
      res = -polylog(m_1, y, z_flat(m_1))!-polylog(m_1,y/z_flat(m_1))
368 369 370
      return
    end if

Luca's avatar
Luca committed
371
    ! need remove trailing zeroes?
372 373 374
    k = size(z_flat)
    kminusj = find_amount_trailing_zeros(z_flat)
    j = k - kminusj
ulrich_y's avatar
ulrich_y committed
375
    if(kminusj > 0) then
ulrich_y's avatar
ulrich_y committed
376
#ifdef DEBUG
Luca's avatar
Luca committed
377
      if(verb >= 30) print*, 'need to remove trailing zeroes'
ulrich_y's avatar
ulrich_y committed
378
#endif
379 380
      allocate(s(j,j))
      s = shuffle_with_zero(z_flat(1:j-1))
ulrich_y's avatar
ulrich_y committed
381
      res = GPL_zero_zi(1,y)*G_flat(z_flat(1:size(z_flat)-1),y)
382
      do i = 1,size(s,1)
Luca's avatar
Luca committed
383
        res = res - G_flat([s(i,:),z_flat(j),zeroes(kminusj-1)], y)
384 385 386 387 388 389
      end do
      res = res / kminusj
      deallocate(s)
      return
    end if

ulrich_y's avatar
ulrich_y committed
390 391 392 393 394 395 396 397 398
    ! is complex?
    if(abs(aimag(y)).gt.zero) then
#ifdef DEBUG
      if(verb >= 30) print*, 'argument is complex, normalising'
#endif
      ! Here we *assume* that y is positive and doesn't mess up the
      ! ieps, which is what GiNaC does (l. 1013)
      do j=1,size(z_flat)
        znorm(j) = inum(z_flat(j)%c/y%c, z_flat(j)%i0)
ulrich_y's avatar
ulrich_y committed
399 400 401
        if (abs(aimag(znorm(j)))>zero) then
          znorm(j)%i0 = int(sign(1._prec, aimag(znorm(j))),1)
        endif
ulrich_y's avatar
ulrich_y committed
402 403 404 405 406
      enddo
      res = G_flat(znorm,inum((1.,0.), y%i0))
      return
    endif

Luca's avatar
Luca committed
407 408
    ! need make convergent?
    if(.not. is_convergent(z_flat,y)) then
ulrich_y's avatar
ulrich_y committed
409
#ifdef DEBUG
Luca's avatar
Luca committed
410
      if(verb >= 10) print*, 'need to make convergent'
ulrich_y's avatar
ulrich_y committed
411
#endif
Luca's avatar
Luca committed
412
      res = make_convergent(z_flat, y)
Luca's avatar
Luca committed
413 414 415
      return
    end if

Luca Naterop's avatar
Luca Naterop committed
416
    ! requires Hoelder convolution?
ulrich_y's avatar
ulrich_y committed
417
    if( any(1.0 <= abs(z_flat%c/y%c) .and. abs(z_flat%c/y%c) <= HoelderCircle) ) then
ulrich_y's avatar
ulrich_y committed
418
      ! Here we *assume* that y is positive and doesn't mess up the
ulrich_y's avatar
ulrich_y committed
419 420 421
      ! ieps, which is what GiNaC does (l. 1013)
      do j=1,size(z_flat)
        znorm(j) = inum(z_flat(j)%c/y%c, z_flat(j)%i0)
ulrich_y's avatar
ulrich_y committed
422 423 424
        if (abs(aimag(znorm(j)))>zero) then
          znorm(j)%i0 = int(sign(1._prec, aimag(znorm(j))),1)
        endif
ulrich_y's avatar
ulrich_y committed
425 426
      enddo
      res = improve_convergence(znorm)
Luca Naterop's avatar
Luca Naterop committed
427 428 429
      return
    end if

Luca's avatar
Luca committed
430
    ! thus it is convergent, and has no trailing zeros
431
    ! -> evaluate in condensed notation which will give series representation
432 433 434 435 436 437 438 439 440 441
    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)
442
    res = G_condensed(m,z,y,size(m))
443 444
    deallocate(m)
    deallocate(z)
445
  END FUNCTION G_flat
446

447 448 449
  FUNCTION G_superflat(g) result(res)
    ! simpler notation for flat evaluation
    complex(kind=prec) :: g(:), res
ulrich_y's avatar
ulrich_y committed
450
    res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g))))
451 452
  END FUNCTION G_superflat

453 454 455 456
  FUNCTION G_real(g) result(res)
    ! simpler notation for flat evaluation
    real(kind=prec) :: g(:)
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
457
    res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g))))
458 459 460 461 462 463
  END FUNCTION G_real

  FUNCTION G_int(g) result(res)
    ! simpler notation for flat evaluation
    integer:: g(:)
    complex(kind=prec) :: res
ulrich_y's avatar
ulrich_y committed
464
    res = G_flat(toinum(real(g(1:size(g)-1),kind=prec)), toinum(real(g(size(g)),kind=prec)))
465 466
  END FUNCTION G_int

467
  RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res)
Luca Naterop's avatar
Luca Naterop committed
468
    ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y)
Luca's avatar
Luca committed
469
    ! assumes zero arguments expressed through the m's
Luca's avatar
Luca committed
470
    
471
    integer :: m(:), k, i
ulrich_y's avatar
ulrich_y committed
472 473
    type(inum) :: z(:), y, z_flat(sum(m))
    complex(kind=prec) :: res, x(k)
474

475 476 477
    ! print*, 'called G_condensed with args'
    ! print*, 'm = ', m
    ! print*, 'z = ', abs(z)
478

479
    ! are all z_i = 0 ? 
480
    if(k == 1 .and. abs(z(1)) < zero) then
481 482
      ! assumes that the zeros at the beginning are passed through m1
      res = GPL_zero_zi(m(1),y)
Luca's avatar
Luca committed
483 484
      return
    end if
485

486 487 488
    ! has trailing zeroes?
    if(abs(z(k)) < zero ) then
      ! we remove them in flat form
489
      z_flat = get_flattened_z(m,z)
490
      res = G_flat(z_flat,y)
491
      return
492 493
    end  if

Luca's avatar
Luca committed
494
    ! need make convergent?
495
    if(.not. all(abs(y) <= abs(z))) then
496 497
      z_flat = get_flattened_z(m,z)
      res = G_flat(z_flat,y)
498
      return
499
    end if
ulrich_y's avatar
ulrich_y committed
500
    ! This is okay because in the MPLs the ieps doesn't matter anymore
ulrich_y's avatar
ulrich_y committed
501
    x(1) = y%c/z(1)%c
ulrich_y's avatar
ulrich_y committed
502
    do i = 2,k
ulrich_y's avatar
ulrich_y committed
503
      x(i) = z(i-1)%c/z(i)%c
Luca Naterop's avatar
Luca Naterop committed
504
    end do
Luca's avatar
Luca committed
505 506 507
    ! print*, 'computed using Li with '
    ! print*, 'm = ', m
    ! print*, 'x = ', x
508
    res = (-1)**k * MPL(m,x)
509
  END FUNCTION G_condensed
Luca's avatar
Luca committed
510

ulrich_y's avatar
ulrich_y committed
511 512 513 514


  FUNCTION G_SUPERFLATN(c0,n)
    integer, intent(in) :: n
ulrich_y's avatar
ulrich_y committed
515
    type(inum), intent(in) :: c0(n)
ulrich_y's avatar
ulrich_y committed
516
    complex(kind=prec) g_superflatn
ulrich_y's avatar
ulrich_y committed
517
    G_superflatn = G_flat(c0(1:n-1), c0(n))
ulrich_y's avatar
ulrich_y committed
518 519
  END FUNCTION

ulrich_y's avatar
ulrich_y committed
520 521 522
  FUNCTION G_FLATr(Z_FLAT,Y)
    real(kind=prec), intent(in) :: z_flat(:), y
    complex(kind=prec) :: g_flatr
ulrich_y's avatar
ulrich_y committed
523
    g_flatr = G_flat(toinum(z_flat), toinum(y))
ulrich_y's avatar
ulrich_y committed
524 525 526 527 528
  END FUNCTION

  FUNCTION G_FLATc(Z_FLAT,Y)
    complex(kind=prec), intent(in) :: z_flat(:), y
    complex(kind=prec) :: g_flatc
ulrich_y's avatar
ulrich_y committed
529
    g_flatc = G_flat(toinum(z_flat), toinum(y))
ulrich_y's avatar
ulrich_y committed
530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567
  END FUNCTION


  FUNCTION G_CONDENSEDs(M,Z,Y)
    implicit none
    type(inum), intent(in) :: z(:), y
    integer m(:)
    complex(kind=prec) g_condenseds
    if (size(m) .ne. size(z)) then
      print*,"G_condesed: weight and args must have the same length",size(m),size(z)
      stop
    endif
    g_condensedS = G_condensed(m, z, y, size(m))
  END FUNCTION
  FUNCTION G_CONDENSEDr(M,Z,Y)
    implicit none
    real(kind=prec), intent(in) :: z(:), y
    integer m(:)
    complex(kind=prec) g_condensedr
    if (size(m) .ne. size(z)) then
      print*,"G_condesed: weight and args must have the same length",size(m),size(z)
      stop
    endif
    g_condensedr = G_condensed(m, toinum(z), inum(y,di0), size(m))
  END FUNCTION
  FUNCTION G_CONDENSEDc(M,Z,Y)
    implicit none
    complex(kind=prec), intent(in) :: z(:), y
    integer m(:)
    complex(kind=prec) g_condensedc
    if (size(m) .ne. size(z)) then
      print*,"G_condesed: weight and args must have the same length",size(m),size(z)
      stop
    endif
    g_condensedc = G_condensed(m, toinum(z), inum(y,di0), size(m))
  END FUNCTION


Luca's avatar
Luca committed
568
END MODULE gpl_module
Luca Naterop's avatar
Luca Naterop committed
569