gpl_module.f90 17.4 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
      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_)
131 132
      return
    end if
Luca's avatar
Luca committed
133
  
134
    a = g(1:size(g)-1)
Luca Naterop's avatar
Luca Naterop committed
135 136 137
    y2 = g(size(g)) 
    m = size(g)  ! actually, size(g)-1+1

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

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

Luca Naterop's avatar
minor  
Luca Naterop committed
196 197 198 199 200
    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
201 202
  END FUNCTION pending_integral

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

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

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

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

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

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

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

ulrich_y's avatar
ulrich_y committed
255
      res = G_flat([izero, a(i+1:size(a))], y2) &
256 257
        + 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
258
        - G_flat([a(i+1)],sr) * G_flat(a(i+1:size(a)), y2)
Luca's avatar
Luca committed
259 260
      return
    end if
Luca's avatar
Luca committed
261

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

ulrich_y's avatar
ulrich_y committed
278
    res = G_flat([a(1:i-1),izero,a(i+1:size(a))],y2) &
Luca's avatar
Luca committed
279 280 281 282
      - 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
283
  END FUNCTION make_convergent
Luca's avatar
Luca committed
284

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

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

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

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

327

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

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

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

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

Luca's avatar
Luca committed
389 390
    ! need make convergent?
    if(.not. is_convergent(z_flat,y)) then
ulrich_y's avatar
ulrich_y committed
391
#ifdef DEBUG
Luca's avatar
Luca committed
392
      if(verb >= 10) print*, 'need to make convergent'
ulrich_y's avatar
ulrich_y committed
393
#endif
Luca's avatar
Luca committed
394
      res = make_convergent(z_flat, y)
Luca's avatar
Luca committed
395 396 397
      return
    end if

Luca Naterop's avatar
Luca Naterop committed
398
    ! requires Hoelder convolution?
ulrich_y's avatar
ulrich_y committed
399
    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
400
      ! Here we *assume* that y is positive and doesn't mess up the
ulrich_y's avatar
ulrich_y committed
401 402 403 404 405
      ! 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
406 407 408
      return
    end if

Luca's avatar
Luca committed
409
    ! thus it is convergent, and has no trailing zeros
410
    ! -> evaluate in condensed notation which will give series representation
411 412 413 414 415 416 417 418 419 420
    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)
421
    res = G_condensed(m,z,y,size(m))
422 423
    deallocate(m)
    deallocate(z)
424
  END FUNCTION G_flat
425

426 427 428
  FUNCTION G_superflat(g) result(res)
    ! simpler notation for flat evaluation
    complex(kind=prec) :: g(:), res
ulrich_y's avatar
ulrich_y committed
429
    res = G_flat(toinum(g(1:size(g)-1)), inum(g(size(g)),di0))
430 431
  END FUNCTION G_superflat

432 433 434 435
  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
436
    res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g))))
437 438 439 440 441 442
  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
443
    res = G_flat(toinum(real(g(1:size(g)-1),kind=prec)), toinum(real(g(size(g)),kind=prec)))
444 445
  END FUNCTION G_int

446
  RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res)
Luca Naterop's avatar
Luca Naterop committed
447
    ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y)
Luca's avatar
Luca committed
448
    ! assumes zero arguments expressed through the m's
Luca's avatar
Luca committed
449
    
450
    integer :: m(:), k, i
ulrich_y's avatar
ulrich_y committed
451 452
    type(inum) :: z(:), y, z_flat(sum(m))
    complex(kind=prec) :: res, x(k)
453

454 455 456
    ! print*, 'called G_condensed with args'
    ! print*, 'm = ', m
    ! print*, 'z = ', abs(z)
457

458
    ! are all z_i = 0 ? 
459
    if(k == 1 .and. abs(z(1)) < zero) then
460 461
      ! assumes that the zeros at the beginning are passed through m1
      res = GPL_zero_zi(m(1),y)
Luca's avatar
Luca committed
462 463
      return
    end if
464

465 466 467
    ! has trailing zeroes?
    if(abs(z(k)) < zero ) then
      ! we remove them in flat form
468
      z_flat = get_flattened_z(m,z)
469
      res = G_flat(z_flat,y)
470
      return
471 472
    end  if

Luca's avatar
Luca committed
473
    ! need make convergent?
474
    if(.not. all(abs(y) <= abs(z))) then
475 476
      z_flat = get_flattened_z(m,z)
      res = G_flat(z_flat,y)
477
      return
478
    end if
ulrich_y's avatar
ulrich_y committed
479
    ! This is okay because in the MPLs the ieps doesn't matter anymore
ulrich_y's avatar
ulrich_y committed
480
    x(1) = y%c/z(1)%c
ulrich_y's avatar
ulrich_y committed
481
    do i = 2,k
ulrich_y's avatar
ulrich_y committed
482
      x(i) = z(i-1)%c/z(i)%c
Luca Naterop's avatar
Luca Naterop committed
483
    end do
Luca's avatar
Luca committed
484 485 486
    ! print*, 'computed using Li with '
    ! print*, 'm = ', m
    ! print*, 'x = ', x
487
    res = (-1)**k * MPL(m,x)
488
  END FUNCTION G_condensed
Luca's avatar
Luca committed
489

ulrich_y's avatar
ulrich_y committed
490 491 492 493


  FUNCTION G_SUPERFLATN(c0,n)
    integer, intent(in) :: n
ulrich_y's avatar
ulrich_y committed
494
    type(inum), intent(in) :: c0(n)
ulrich_y's avatar
ulrich_y committed
495
    complex(kind=prec) g_superflatn
ulrich_y's avatar
ulrich_y committed
496
    G_superflatn = G_flat(c0(1:n-1), c0(n))
ulrich_y's avatar
ulrich_y committed
497 498
  END FUNCTION

ulrich_y's avatar
ulrich_y committed
499 500 501
  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
502
    g_flatr = G_flat(toinum(z_flat), toinum(y))
ulrich_y's avatar
ulrich_y committed
503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 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
  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
547
END MODULE gpl_module
Luca Naterop's avatar
Luca Naterop committed
548