gpl_module.f90 17.9 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

55
  RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) 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 59
    type(inum) :: a(:), y2, s(m), p(:)
    complex(kind=prec) :: res
60
    type(inum) :: alpha(product((/(i,i=1,size(a)+size(s))/))/  &
61 62 63
        (product((/(i,i=1,size(a))/))*product((/(i,i=1,size(s))/))), & 
        size(a) + size(s))

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

88
  RECURSIVE FUNCTION pending_integral(p,i,g) result(res)
89
    ! evaluates a pending integral by reducing it to simpler ones and g functions
ulrich_y's avatar
ulrich_y committed
90 91 92
    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
93
    integer :: i, m
94
    res = 0
95

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

    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
105
#ifdef DEBUG
106
      if(verb >= 30) print*, 'only integrals in front, we get G-function'
ulrich_y's avatar
ulrich_y committed
107
#endif
108 109 110
      res = G_flat(b,y1)
      return
    end if
Luca's avatar
Luca committed
111

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

122
    ! if depth one and m = 1 use (23) 
123
    if(size(g) == 1) then
ulrich_y's avatar
ulrich_y committed
124
#ifdef DEBUG
125
      if(verb >= 30) print*, 'case depth one and m = 1'
ulrich_y's avatar
ulrich_y committed
126
#endif
ulrich_y's avatar
ulrich_y committed
127 128
      !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
129 130 131
      if(abs(aimag(p(1))).gt.zero) then
        p(1)%i0 = int(sign(1._prec, aimag(p(1))),1)
      endif
ulrich_y's avatar
ulrich_y committed
132 133
      res = pending_integral(p,2,[inum( g(1)%c,-g(1)%i0 )  ]) - pending_integral(p,2,[izero]) &
        + G_flat(p(2:size(p)), p(1)) * (log(g(1)%c) + p(1)%i0 * pi * i_)
134 135
      return
    end if
Luca's avatar
Luca committed
136
  
137
    a = g(1:size(g)-1)
Luca Naterop's avatar
Luca Naterop committed
138 139 140
    y2 = g(size(g)) 
    m = size(g)  ! actually, size(g)-1+1

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

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

Luca Naterop's avatar
minor  
Luca Naterop committed
199 200 201 202 203
    res =  +pending_integral(p,1,zeroes(0)) * 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]) & 
      + pending_integral([p,a(i-1)],1,zeroes(0)) * G_flat(a,y2) & 
      + pending_integral([p,a(i)], i, [a(1:i-1), a(i+1:size(a)),y2])  & 
      - pending_integral([p,a(i)],1,zeroes(0)) * G_flat(a,y2) 
Luca's avatar
Luca committed
204 205
  END FUNCTION pending_integral

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

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

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

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

Luca Naterop's avatar
Luca Naterop committed
231
    if(i == size(a)) then
232
      ! sr at the end, thus shuffle
ulrich_y's avatar
ulrich_y committed
233
#ifdef DEBUG
Luca Naterop's avatar
Luca Naterop committed
234
      if(verb >= 30) print*, 'sr at the end'
ulrich_y's avatar
ulrich_y committed
235
#endif
Luca Naterop's avatar
Luca Naterop committed
236 237 238 239
      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
240

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

ulrich_y's avatar
ulrich_y committed
244
#ifdef DEBUG
Luca's avatar
Luca committed
245 246 247
      if(verb >= 30) then 
        print*, '--------------------------------------------------'
        print*, 'sr at beginning, map to: '
ulrich_y's avatar
ulrich_y committed
248
        call print_G([izero, a(i+1:size(a))], y2)
Luca's avatar
Luca committed
249 250
        call print_G([y2], sr)
        call print_G(a(i+1:size(a)), y2)
251
        print*, 'PI with p=',real([sr, a(i+1)]),'i=',i,'g =', real([a(i+2:size(a)), y2])
Luca's avatar
Luca committed
252 253 254 255
        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
256
#endif
Luca's avatar
Luca committed
257

ulrich_y's avatar
ulrich_y committed
258
      res = G_flat([izero, a(i+1:size(a))], y2) &
259 260
        + G_flat([y2], sr) * G_flat(a(i+1:size(a)), y2) &
        + pending_integral([sr, a(i+1)], i, [a(i+2:size(a)), y2]) &
ulrich_y's avatar
ulrich_y committed
261
        - G_flat([a(i+1)],sr) * G_flat(a(i+1:size(a)), y2)
Luca's avatar
Luca committed
262 263
      return
    end if
Luca's avatar
Luca committed
264

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

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

ulrich_y's avatar
ulrich_y committed
288
  RECURSIVE FUNCTION improve_convergence(z) result(res)
Luca Naterop's avatar
Luca Naterop committed
289
    ! improves the convergence by applying the Hoelder convolution to G(z1,...zk,1)
ulrich_y's avatar
ulrich_y committed
290 291
    type(inum) :: z(:),oneminusz(size(z))
    complex(kind=prec) :: res
Luca Naterop's avatar
Luca Naterop committed
292 293
    complex(kind=prec), parameter :: p = 2.0
    integer :: k, j
ulrich_y's avatar
ulrich_y committed
294
#ifdef DEBUG
Luca Naterop's avatar
Luca Naterop committed
295
    if(verb >= 30) print*, 'requires Hoelder convolution'
ulrich_y's avatar
ulrich_y committed
296
#endif
ulrich_y's avatar
ulrich_y committed
297 298 299 300
    ! 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
301
    do j=1,size(z)
ulrich_y's avatar
ulrich_y committed
302 303 304 305 306
      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
307
    enddo
Luca Naterop's avatar
Luca Naterop committed
308 309
    k = size(z)

ulrich_y's avatar
ulrich_y committed
310 311
    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
312
    do j = 1,k-1
ulrich_y's avatar
ulrich_y committed
313
      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
314 315 316
    end do
  END FUNCTION improve_convergence

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

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

330

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

338
    ! add small imaginary part if not there
Luca Naterop's avatar
Luca Naterop committed
339 340 341 342
    ! 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
343

344
    ! is just a logarithm? 
Luca's avatar
Luca committed
345
    if(all(abs(z_flat) < zero)) then
ulrich_y's avatar
ulrich_y committed
346
#ifdef DEBUG
347
      if(verb >= 70) print*, 'all z are zero'
ulrich_y's avatar
ulrich_y committed
348
#endif
ulrich_y's avatar
ulrich_y committed
349
      res = gpl_zero_zi(size(z_flat),y)
Luca's avatar
Luca committed
350 351
      return
    end if
352
    if(size(z_flat) == 1) then
ulrich_y's avatar
ulrich_y committed
353
#ifdef DEBUG
354
      if(verb >= 70) print*, 'is just a logarithm'
ulrich_y's avatar
ulrich_y committed
355
#endif
ulrich_y's avatar
ulrich_y committed
356
      res = plog1(y,z_flat(1))  ! log((z_flat(1) - y)/z_flat(1))
357 358 359 360 361 362 363 364 365
      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
366
#ifdef DEBUG
Luca's avatar
Luca committed
367
      if(verb >= 70) print*, 'is just a polylog'
ulrich_y's avatar
ulrich_y committed
368
#endif
ulrich_y's avatar
ulrich_y committed
369
      res = -polylog(m_1, y, z_flat(m_1))!-polylog(m_1,y/z_flat(m_1))
370 371 372
      return
    end if

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

ulrich_y's avatar
ulrich_y committed
392 393 394 395 396 397 398 399 400 401 402 403 404 405
    ! 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)
      enddo
      res = G_flat(znorm,inum((1.,0.), y%i0))
      return
    endif

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

Luca Naterop's avatar
Luca Naterop committed
415
    ! requires Hoelder convolution?
ulrich_y's avatar
ulrich_y committed
416
    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
417
      ! Here we *assume* that y is positive and doesn't mess up the
ulrich_y's avatar
ulrich_y committed
418 419 420 421 422
      ! 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)
      enddo
      res = improve_convergence(znorm)
Luca Naterop's avatar
Luca Naterop committed
423 424 425
      return
    end if

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

443 444 445
  FUNCTION G_superflat(g) result(res)
    ! simpler notation for flat evaluation
    complex(kind=prec) :: g(:), res
ulrich_y's avatar
ulrich_y committed
446
    res = G_flat(toinum(g(1:size(g)-1)), inum(g(size(g)),di0))
447 448
  END FUNCTION G_superflat

449 450 451 452
  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
453
    res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g))))
454 455 456 457 458 459
  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
460
    res = G_flat(toinum(real(g(1:size(g)-1),kind=prec)), toinum(real(g(size(g)),kind=prec)))
461 462
  END FUNCTION G_int

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

471 472 473
    ! print*, 'called G_condensed with args'
    ! print*, 'm = ', m
    ! print*, 'z = ', abs(z)
474

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

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

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

ulrich_y's avatar
ulrich_y committed
507 508 509 510


  FUNCTION G_SUPERFLATN(c0,n)
    integer, intent(in) :: n
ulrich_y's avatar
ulrich_y committed
511
    type(inum), intent(in) :: c0(n)
ulrich_y's avatar
ulrich_y committed
512
    complex(kind=prec) g_superflatn
ulrich_y's avatar
ulrich_y committed
513
    G_superflatn = G_flat(c0(1:n-1), c0(n))
ulrich_y's avatar
ulrich_y committed
514 515
  END FUNCTION

ulrich_y's avatar
ulrich_y committed
516 517 518
  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
519
    g_flatr = G_flat(toinum(z_flat), toinum(y))
ulrich_y's avatar
ulrich_y committed
520 521 522 523 524 525 526 527 528 529 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
  END FUNCTION

  FUNCTION G_FLATc(Z_FLAT,Y)
    complex(kind=prec), intent(in) :: z_flat(:), y
    complex(kind=prec) :: g_flatc
    g_flatc = G_flat(toinum(z_flat), inum(y,di0))
  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
564
END MODULE gpl_module
Luca Naterop's avatar
Luca Naterop committed
565