gpl_module.f90 18.2 KB
 Luca Naterop committed Apr 04, 2019 1 `````` `````` Luca committed Apr 04, 2019 2 ``````MODULE gpl_module `````` Luca committed May 08, 2019 3 `````` use globals `````` ulrich_y committed Jul 08, 2019 4 `````` use ieps `````` ulrich_y committed Aug 20, 2019 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 committed Apr 04, 2019 9 10 `````` implicit none `````` ulrich_y committed Jul 10, 2019 11 `````` INTERFACE G `````` ulrich_y committed Jul 19, 2019 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 committed Jul 10, 2019 15 `````` END INTERFACE G `````` Luca committed Apr 04, 2019 16 17 ``````CONTAINS `````` Luca committed Apr 09, 2019 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 committed Jul 08, 2019 21 22 `````` type(inum) :: y complex(kind=prec) :: GPL_zero_zi `````` ulrich_y committed Jul 09, 2019 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 committed Aug 22, 2019 27 `````` GPL_zero_zi = 1.0_prec/factorial(l) * (log(-real(y))+i_*(y%i0*pi)) ** l `````` ulrich_y committed Jul 09, 2019 28 29 30 31 `````` endif else GPL_zero_zi = 1.0_prec/factorial(l) * log(y%c) ** l endif `````` Luca committed Apr 09, 2019 32 33 `````` END FUNCTION GPL_zero_zi `````` Luca committed May 06, 2019 34 35 `````` FUNCTION is_convergent(z,y) ! returns true if G(z,y) convergent, otherwise false `````` Luca committed May 15, 2019 36 `````` ! can be used in either notation (flat or condensed) `````` ulrich_y committed Jul 08, 2019 37 `````` type(inum) :: z(:), y `````` Luca committed May 06, 2019 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 committed May 15, 2019 48 `````` SUBROUTINE print_G(z_flat, y) `````` ulrich_y committed Jul 08, 2019 49 50 `````` type(inum) :: z_flat(:) type(inum) , optional :: y `````` Luca Naterop committed Jul 08, 2019 51 `````` if(present(y)) print*, 'G(', real(z_flat), real(y), ')' `````` Luca committed May 23, 2019 52 `````` if(.not. present(y)) print*, 'G(', abs(z_flat), ')' `````` Luca committed May 15, 2019 53 54 `````` END SUBROUTINE print_G `````` ulrich_y committed Sep 27, 2019 55 `````` RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p, srs) result(res) `````` Luca Naterop committed Jul 05, 2019 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 committed Sep 27, 2019 58 `````` integer(1) :: srs `````` ulrich_y committed Jul 08, 2019 59 60 `````` type(inum) :: a(:), y2, s(m), p(:) complex(kind=prec) :: res `````` ulrich_y committed Jul 08, 2019 61 `````` type(inum) :: alpha(product((/(i,i=1,size(a)+size(s))/))/ & `````` Luca Naterop committed Jul 05, 2019 62 63 64 `````` (product((/(i,i=1,size(a))/))*product((/(i,i=1,size(s))/))), & size(a) + size(s)) `````` ulrich_y committed Jul 08, 2019 65 `````` s = [zeroes(m-1),marker] `````` Luca Naterop committed Jul 05, 2019 66 `````` alpha = shuffle_product(a,s) `````` ulrich_y committed Jul 10, 2019 67 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 68 69 70 `````` if(verb >= 50) then print*, 'mapping to ' call print_G(a,y2) `````` Luca Naterop committed Jul 08, 2019 71 `````` print*, 'PI with p=',real(p),'i=',m,'g =',real([zeroes(m-1),y2]) `````` Luca Naterop committed Jul 05, 2019 72 `````` end if `````` ulrich_y committed Jul 10, 2019 73 ``````#endif `````` ulrich_y committed Sep 27, 2019 74 `````` res = G_flat(a,y2)*pending_integral(p,m,[zeroes(m-1),y2], srs) `````` ulrich_y committed Jul 10, 2019 75 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 76 `````` if(verb >= 50) print*, 'also mapping to' `````` ulrich_y committed Jul 10, 2019 77 ``````#endif `````` Luca Naterop committed Jul 05, 2019 78 79 `````` do j = 2,size(alpha, 1) ! find location of s_r `````` ulrich_y committed Jul 08, 2019 80 `````` n = find_marker(alpha(j,:)) `````` ulrich_y committed Jul 10, 2019 81 ``````#ifdef DEBUG `````` ulrich_y committed Jul 08, 2019 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 committed Jul 10, 2019 84 ``````#endif `````` ulrich_y committed Sep 27, 2019 85 `````` res = res - pending_integral(p, n, [alpha(j,1:n-1),alpha(j,n+1:size(alpha,2)),y2], srs) `````` Luca Naterop committed Jul 05, 2019 86 87 88 `````` end do END FUNCTION remove_sr_from_last_place_in_PI `````` ulrich_y committed Sep 27, 2019 89 `````` RECURSIVE FUNCTION pending_integral(p,i,g,srs) result(res) `````` Luca committed May 21, 2019 90 `````` ! evaluates a pending integral by reducing it to simpler ones and g functions `````` ulrich_y committed Jul 08, 2019 91 92 93 `````` complex(kind=prec) :: res type(inum) :: p(:), g(:) type(inum) :: y1, y2, b(size(p)-1), a(size(g)-1) `````` Luca committed Jun 04, 2019 94 `````` integer :: i, m `````` ulrich_y committed Sep 27, 2019 95 `````` integer(1):: srs `````` Luca committed May 21, 2019 96 `````` res = 0 `````` Luca committed May 23, 2019 97 `````` `````` ulrich_y committed Jul 10, 2019 98 ``````#ifdef DEBUG `````` ulrich_y committed Jul 08, 2019 99 `````` if(verb >= 30) print*, 'evaluating PI with p=',real(p),'i=',real(i),'g =',real(g) `````` ulrich_y committed Jul 10, 2019 100 ``````#endif `````` Luca Naterop committed Jul 05, 2019 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 committed Jul 10, 2019 107 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 108 `````` if(verb >= 30) print*, 'only integrals in front, we get G-function' `````` ulrich_y committed Jul 10, 2019 109 ``````#endif `````` Luca Naterop committed Jul 05, 2019 110 111 112 `````` res = G_flat(b,y1) return end if `````` Luca committed Jun 04, 2019 113 `````` `````` Luca committed May 23, 2019 114 `````` ! if integration variable at end -> we gat a G function `````` Luca committed May 21, 2019 115 `````` if(i == size(g)+1) then `````` ulrich_y committed Jul 10, 2019 116 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 117 `````` if(verb >= 30) print*, 'is just a G-function' `````` ulrich_y committed Jul 10, 2019 118 ``````#endif `````` Luca committed May 21, 2019 119 120 121 `````` res = G_flat([p(2:size(p)),g], p(1)) return end if `````` Luca committed Jun 04, 2019 122 `````` `````` Luca Naterop committed Jul 05, 2019 123 `````` `````` Luca Naterop committed Aug 28, 2019 124 `````` ! if depth one and m = 1 use (23) `````` Luca committed May 21, 2019 125 `````` if(size(g) == 1) then `````` ulrich_y committed Jul 10, 2019 126 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 127 `````` if(verb >= 30) print*, 'case depth one and m = 1' `````` ulrich_y committed Jul 10, 2019 128 ``````#endif `````` ulrich_y committed Jul 08, 2019 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 committed Sep 27, 2019 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_) `````` Luca committed May 21, 2019 133 134 `````` return end if `````` Luca committed Jun 04, 2019 135 `````` `````` ulrich_y committed Jul 09, 2019 136 `````` a = g(1:size(g)-1) `````` Luca Naterop committed Jul 05, 2019 137 138 139 `````` y2 = g(size(g)) m = size(g) ! actually, size(g)-1+1 `````` Luca Naterop committed Aug 28, 2019 140 `````` ! if depth one and m > 1 use (25) `````` Luca Naterop committed Jul 05, 2019 141 `````` if(all( abs( g(1:size(g)-1) ) < zero)) then `````` ulrich_y committed Jul 10, 2019 142 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 143 144 145 146 `````` if(verb >= 30) print*, 'case depth one and m > 1' if(verb >= 50) then print*, 'map to' print*, 'zeta(',m,')' `````` ulrich_y committed Jul 08, 2019 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 committed Jul 08, 2019 150 `````` print*, 'PI with p=',tocmplx([p,izero]),'i=',m-1,'g =',tocmplx(zeroes(0)) `````` Luca Naterop committed Jul 05, 2019 151 `````` end if `````` ulrich_y committed Jul 10, 2019 152 ``````#endif `````` ulrich_y committed Sep 27, 2019 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 committed Jun 04, 2019 156 157 158 `````` return end if `````` Luca Naterop committed Aug 28, 2019 159 `````` ! case of higher depth, s_r at beginning, use (21) `````` Luca committed May 23, 2019 160 `````` if(i == 1) then `````` ulrich_y committed Jul 10, 2019 161 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 162 `````` if(verb >= 30) print*, 'case higher depth, sr at beginning' `````` Luca committed Jun 04, 2019 163 `````` `````` Luca Naterop committed Jul 05, 2019 164 `````` if(verb >= 50) then `````` Luca Naterop committed Aug 28, 2019 165 `````` print*, 'map to (using 21)' `````` ulrich_y committed Jul 08, 2019 166 `````` print*, 'PI with p=',real(p),'i=',0,'g =',tocmplx(zeroes(0) ) `````` ulrich_y committed Jul 08, 2019 167 `````` call print_G([izero,a],y2) `````` ulrich_y committed Jul 08, 2019 168 `````` print*, 'PI with p=',real([p,y2]),'i=',0,'g =',tocmplx(zeroes(0) ) `````` Luca Naterop committed Jul 05, 2019 169 `````` call print_G(a,y2) `````` Luca Naterop committed Jul 08, 2019 170 `````` print*, 'PI with p=',real([p,a(1)]),'i=',1,'g =',real([a(2:size(a)),y2]) `````` ulrich_y committed Jul 08, 2019 171 `````` print*, 'PI with p=',real([p,a(1)]),'i=',0,'g =',tocmplx(zeroes(0)) `````` Luca Naterop committed Jul 05, 2019 172 173 `````` call print_G(a,y2) end if `````` ulrich_y committed Jul 10, 2019 174 ``````#endif `````` Luca Naterop committed Jul 05, 2019 175 `````` `````` ulrich_y committed Sep 27, 2019 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) `````` Luca committed May 23, 2019 180 181 `````` return end if `````` Luca committed Jun 04, 2019 182 `````` `````` Luca Naterop committed Aug 28, 2019 183 `````` ! case higher depth, s_r at the end, use shuffle algebra to remove from last place `````` Luca committed Jun 04, 2019 184 `````` if(i == size(g)) then `````` ulrich_y committed Jul 10, 2019 185 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 186 `````` if(verb >= 30) print*, 's_r at the end under PI, need to shuffle' `````` ulrich_y committed Jul 10, 2019 187 ``````#endif `````` Luca Naterop committed Jul 05, 2019 188 `````` m = find_amount_trailing_zeros(a) + 1 `````` ulrich_y committed Sep 27, 2019 189 `````` res = remove_sr_from_last_place_in_PI(a(1:size(a)-(m-1)), y2, m, p, srs) `````` Luca committed Jun 04, 2019 190 191 192 `````` return end if `````` Luca Naterop committed Aug 28, 2019 193 `````` ! case higher depth, s_r in middle, use (22) `````` ulrich_y committed Jul 10, 2019 194 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 195 `````` if(verb >= 30) print*, 's_r in the middle under PI' `````` ulrich_y committed Jul 10, 2019 196 ``````#endif `````` Luca Naterop committed Jul 07, 2019 197 `````` `````` ulrich_y committed Sep 27, 2019 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 committed May 15, 2019 203 204 `````` END FUNCTION pending_integral `````` Luca Naterop committed Jul 08, 2019 205 `````` RECURSIVE FUNCTION remove_sr_from_last_place_in_G(a,y2,m,sr) result(res) `````` ulrich_y committed Jul 08, 2019 206 207 `````` type(inum) :: a(:), sr, y2 complex(kind=prec) :: res `````` Luca committed May 21, 2019 208 `````` integer :: m,i,j `````` ulrich_y committed Jul 08, 2019 209 `````` type(inum) :: alpha(product((/(i,i=1,size(a)+m)/))/ & `````` Luca committed May 21, 2019 210 211 `````` (product((/(i,i=1,size(a))/))*product((/(i,i=1,m)/))), & size(a) + m) `````` Luca committed Jun 04, 2019 212 213 `````` alpha = shuffle_product(a,[zeroes(m-1),sr]) res = G_flat(a,y2)*G_flat([zeroes(m-1),sr],y2) `````` Luca committed May 21, 2019 214 215 216 `````` do j = 2,size(alpha,1) res = res - G_flat(alpha(j,:),y2) end do `````` Luca Naterop committed Jul 05, 2019 217 `````` END FUNCTION remove_sr_from_last_place_in_G `````` Luca committed May 21, 2019 218 `````` `````` Luca Naterop committed Jul 08, 2019 219 `````` RECURSIVE FUNCTION make_convergent(a,y2) result(res) `````` Luca Naterop committed Aug 28, 2019 220 `````` ! goes from G-functions to pending integrals and simpler expressions `````` Luca committed Jun 04, 2019 221 `````` `````` ulrich_y committed Jul 08, 2019 222 223 `````` type(inum) :: a(:), y2, sr complex(kind=prec) :: res `````` Luca committed May 21, 2019 224 `````` integer :: i, mminus1 `````` Luca committed May 20, 2019 225 `````` `````` Luca committed May 15, 2019 226 `````` res = 0 `````` Luca Naterop committed May 16, 2019 227 228 `````` i = min_index(abs(a)) sr = a(i) `````` Luca committed Jun 04, 2019 229 `````` `````` Luca Naterop committed Jul 07, 2019 230 `````` if(i == size(a)) then `````` Luca Naterop committed Aug 28, 2019 231 `````` ! sr at the end, thus shuffle `````` ulrich_y committed Jul 10, 2019 232 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 233 `````` if(verb >= 30) print*, 'sr at the end' `````` ulrich_y committed Jul 10, 2019 234 ``````#endif `````` Luca Naterop committed Jul 07, 2019 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 committed Jun 04, 2019 239 `````` `````` Luca Naterop committed May 16, 2019 240 `````` if(i == 1) then `````` Luca Naterop committed Aug 28, 2019 241 `````` !s_r at beginning, thus use (18) `````` Luca committed Jun 04, 2019 242 `````` `````` ulrich_y committed Jul 10, 2019 243 ``````#ifdef DEBUG `````` Luca committed Jun 04, 2019 244 245 246 `````` if(verb >= 30) then print*, '--------------------------------------------------' print*, 'sr at beginning, map to: ' `````` ulrich_y committed Jul 08, 2019 247 `````` call print_G([izero, a(i+1:size(a))], y2) `````` Luca committed Jun 04, 2019 248 249 `````` call print_G([y2], sr) call print_G(a(i+1:size(a)), y2) `````` Luca Naterop committed Jul 08, 2019 250 `````` print*, 'PI with p=',real([sr, a(i+1)]),'i=',i,'g =', real([a(i+2:size(a)), y2]) `````` Luca committed Jun 04, 2019 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 committed Jul 10, 2019 255 ``````#endif `````` ulrich_y committed Jul 08, 2019 256 `````` res = G_flat([izero, a(i+1:size(a))], y2) & `````` Luca Naterop committed May 16, 2019 257 `````` + G_flat([y2], sr) * G_flat(a(i+1:size(a)), y2) & `````` ulrich_y committed Sep 27, 2019 258 `````` + pending_integral([sr, a(i+1)], i, [a(i+2:size(a)), y2],sr%i0) & `````` ulrich_y committed Jul 08, 2019 259 `````` - G_flat([a(i+1)],sr) * G_flat(a(i+1:size(a)), y2) `````` Luca committed May 15, 2019 260 261 `````` return end if `````` Luca committed May 21, 2019 262 `````` `````` Luca Naterop committed Aug 28, 2019 263 `````` ! so s_r in middle, use (19) `````` ulrich_y committed Jul 10, 2019 264 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 265 266 267 `````` if(verb >= 30) then print*, '--------------------------------------------------' print*, 'sr in the middle, map to: ' `````` ulrich_y committed Jul 08, 2019 268 `````` call print_G([a(1:i-1),izero,a(i+1:size(a))],y2) `````` Luca Naterop committed Jul 08, 2019 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]) `````` Luca committed May 23, 2019 270 271 `````` call print_G([a(i-1)],sr) call print_G([a(1:i-1),a(i+1:size(a))],y2) `````` Luca Naterop committed Jul 08, 2019 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]) `````` Luca committed May 23, 2019 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 committed Jul 10, 2019 277 ``````#endif `````` Luca committed May 23, 2019 278 `````` `````` ulrich_y committed Jul 08, 2019 279 `````` res = G_flat([a(1:i-1),izero,a(i+1:size(a))],y2) & `````` ulrich_y committed Sep 27, 2019 280 `````` - pending_integral([sr,a(i-1)], i-1, [a(1:i-2),a(i+1:size(a)),y2],sr%i0) & `````` Luca committed May 21, 2019 281 `````` + G_flat([a(i-1)],sr) * G_flat([a(1:i-1),a(i+1:size(a))],y2) & `````` ulrich_y committed Sep 27, 2019 282 `````` + pending_integral([sr,a(i+1)], i, [a(1:i-1),a(i+2:size(a)),y2],sr%i0) & `````` Luca committed May 21, 2019 283 `````` - G_flat([a(i+1)],sr) * G_flat([a(1:i-1),a(i+1:size(a))],y2) `````` Luca committed Jun 04, 2019 284 `````` END FUNCTION make_convergent `````` Luca committed May 15, 2019 285 `````` `````` ulrich_y committed Jul 08, 2019 286 `````` RECURSIVE FUNCTION improve_convergence(z) result(res) `````` Luca Naterop committed Jul 07, 2019 287 `````` ! improves the convergence by applying the Hoelder convolution to G(z1,...zk,1) `````` ulrich_y committed Jul 08, 2019 288 289 `````` type(inum) :: z(:),oneminusz(size(z)) complex(kind=prec) :: res `````` Luca Naterop committed Jul 07, 2019 290 291 `````` complex(kind=prec), parameter :: p = 2.0 integer :: k, j `````` ulrich_y committed Jul 10, 2019 292 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 293 `````` if(verb >= 30) print*, 'requires Hoelder convolution' `````` ulrich_y committed Jul 10, 2019 294 ``````#endif `````` ulrich_y committed Jul 09, 2019 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 committed Jul 09, 2019 299 `````` do j=1,size(z) `````` ulrich_y committed Jul 09, 2019 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 committed Jul 09, 2019 305 `````` enddo `````` Luca Naterop committed Jul 07, 2019 306 307 `````` k = size(z) `````` ulrich_y committed Jul 09, 2019 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 committed Jul 07, 2019 310 `````` do j = 1,k-1 `````` ulrich_y committed Jul 09, 2019 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 committed Jul 07, 2019 312 313 314 `````` end do END FUNCTION improve_convergence `````` 315 `````` RECURSIVE FUNCTION G_flat(z_flat,y) result(res) `````` Luca Naterop committed Apr 29, 2019 316 `````` ! Calls G function with flat arguments, that is, zeroes not passed through the m's. `````` ulrich_y committed Jul 09, 2019 317 `````` type(inum) :: z_flat(:), y, znorm(size(z_flat)) `````` ulrich_y committed Jul 08, 2019 318 319 `````` complex(kind=prec) :: res type(inum), allocatable :: s(:,:), z(:) `````` Luca committed May 20, 2019 320 `````` integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i, m_1 `````` Luca committed Apr 25, 2019 321 `````` integer, allocatable :: m(:) `````` Luca committed May 20, 2019 322 `````` logical :: is_depth_one `````` 323 `````` `````` ulrich_y committed Jul 10, 2019 324 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 325 `````` if(verb >= 50) call print_G(z_flat,y) `````` ulrich_y committed Jul 10, 2019 326 ``````#endif `````` Luca committed May 15, 2019 327 `````` `````` Luca Naterop committed Jul 08, 2019 328 `````` `````` Luca Naterop committed Jul 08, 2019 329 `````` if(size(z_flat) == 1) then `````` ulrich_y committed Jul 09, 2019 330 `````` if( abs(z_flat(1)%c - y%c) <= zero ) then `````` Luca Naterop committed Jul 08, 2019 331 332 333 334 335 `````` res = 0 return end if end if `````` Luca Naterop committed Jul 05, 2019 336 `````` ! add small imaginary part if not there `````` Luca Naterop committed Jul 05, 2019 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 committed May 21, 2019 341 `````` `````` Luca committed May 21, 2019 342 `````` ! is just a logarithm? `````` Luca committed May 21, 2019 343 `````` if(all(abs(z_flat) < zero)) then `````` ulrich_y committed Jul 10, 2019 344 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 345 `````` if(verb >= 70) print*, 'all z are zero' `````` ulrich_y committed Jul 10, 2019 346 ``````#endif `````` ulrich_y committed Jul 09, 2019 347 `````` res = gpl_zero_zi(size(z_flat),y) `````` Luca committed May 21, 2019 348 349 `````` return end if `````` Luca committed May 21, 2019 350 `````` if(size(z_flat) == 1) then `````` ulrich_y committed Jul 10, 2019 351 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 352 `````` if(verb >= 70) print*, 'is just a logarithm' `````` ulrich_y committed Jul 10, 2019 353 ``````#endif `````` ulrich_y committed Jul 09, 2019 354 `````` res = plog1(y,z_flat(1)) ! log((z_flat(1) - y)/z_flat(1)) `````` Luca committed May 21, 2019 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 committed Jul 10, 2019 364 ``````#ifdef DEBUG `````` Luca committed May 22, 2019 365 `````` if(verb >= 70) print*, 'is just a polylog' `````` ulrich_y committed Jul 10, 2019 366 ``````#endif `````` ulrich_y committed Jul 09, 2019 367 `````` res = -polylog(m_1, y, z_flat(m_1))!-polylog(m_1,y/z_flat(m_1)) `````` Luca committed May 21, 2019 368 369 370 `````` return end if `````` Luca committed May 07, 2019 371 `````` ! need remove trailing zeroes? `````` 372 373 374 `````` k = size(z_flat) kminusj = find_amount_trailing_zeros(z_flat) j = k - kminusj `````` ulrich_y committed Aug 19, 2019 375 `````` if(kminusj > 0) then `````` ulrich_y committed Jul 10, 2019 376 ``````#ifdef DEBUG `````` Luca committed May 22, 2019 377 `````` if(verb >= 30) print*, 'need to remove trailing zeroes' `````` ulrich_y committed Jul 10, 2019 378 ``````#endif `````` 379 380 `````` allocate(s(j,j)) s = shuffle_with_zero(z_flat(1:j-1)) `````` ulrich_y committed Jul 09, 2019 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 committed Jun 04, 2019 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 committed Sep 23, 2019 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 committed Sep 27, 2019 399 400 401 `````` if (abs(aimag(znorm(j)))>zero) then znorm(j)%i0 = int(sign(1._prec, aimag(znorm(j))),1) endif `````` ulrich_y committed Sep 23, 2019 402 403 404 405 406 `````` enddo res = G_flat(znorm,inum((1.,0.), y%i0)) return endif `````` Luca committed May 21, 2019 407 408 `````` ! need make convergent? if(.not. is_convergent(z_flat,y)) then `````` ulrich_y committed Jul 10, 2019 409 ``````#ifdef DEBUG `````` Luca committed May 22, 2019 410 `````` if(verb >= 10) print*, 'need to make convergent' `````` ulrich_y committed Jul 10, 2019 411 ``````#endif `````` Luca committed Jun 04, 2019 412 `````` res = make_convergent(z_flat, y) `````` Luca committed May 21, 2019 413 414 415 `````` return end if `````` Luca Naterop committed Jul 07, 2019 416 `````` ! requires Hoelder convolution? `````` ulrich_y committed Jul 10, 2019 417 `````` if( any(1.0 <= abs(z_flat%c/y%c) .and. abs(z_flat%c/y%c) <= HoelderCircle) ) then `````` ulrich_y committed Aug 28, 2019 418 `````` ! Here we *assume* that y is positive and doesn't mess up the `````` ulrich_y committed Jul 09, 2019 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 committed Sep 27, 2019 422 423 424 `````` if (abs(aimag(znorm(j)))>zero) then znorm(j)%i0 = int(sign(1._prec, aimag(znorm(j))),1) endif `````` ulrich_y committed Jul 09, 2019 425 426 `````` enddo res = improve_convergence(znorm) `````` Luca Naterop committed Jul 07, 2019 427 428 429 `````` return end if `````` Luca committed May 15, 2019 430 `````` ! thus it is convergent, and has no trailing zeros `````` Luca Naterop committed Jul 10, 2019 431 `````` ! -> evaluate in condensed notation which will give series representation `````` Luca committed Apr 25, 2019 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)) `````` Luca committed Apr 25, 2019 443 444 `````` deallocate(m) deallocate(z) `````` Luca Naterop committed May 16, 2019 445 `````` END FUNCTION G_flat `````` Luca committed Apr 25, 2019 446 `````` `````` Luca committed May 22, 2019 447 448 449 `````` FUNCTION G_superflat(g) result(res) ! simpler notation for flat evaluation complex(kind=prec) :: g(:), res `````` ulrich_y committed Sep 27, 2019 450 `````` res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g)))) `````` Luca committed May 22, 2019 451 452 `````` END FUNCTION G_superflat `````` Luca committed May 23, 2019 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 committed Aug 22, 2019 457 `````` res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g)))) `````` Luca committed May 23, 2019 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 committed Aug 22, 2019 464 `````` res = G_flat(toinum(real(g(1:size(g)-1),kind=prec)), toinum(real(g(size(g)),kind=prec))) `````` Luca committed May 23, 2019 465 466 `````` END FUNCTION G_int `````` 467 `````` RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res) `````` Luca Naterop committed Apr 04, 2019 468 `````` ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) `````` Luca committed Apr 09, 2019 469 `````` ! assumes zero arguments expressed through the m's `````` Luca committed Apr 23, 2019 470 `````` `````` 471 `````` integer :: m(:), k, i `````` ulrich_y committed Jul 09, 2019 472 473 `````` type(inum) :: z(:), y, z_flat(sum(m)) complex(kind=prec) :: res, x(k) `````` Luca committed Apr 25, 2019 474 `````` `````` Luca committed May 08, 2019 475 476 477 `````` ! print*, 'called G_condensed with args' ! print*, 'm = ', m ! print*, 'z = ', abs(z) `````` 478 `````` `````` Luca committed Apr 09, 2019 479 `````` ! are all z_i = 0 ? `````` 480 `````` if(k == 1 .and. abs(z(1)) < zero) then `````` Luca committed Apr 25, 2019 481 482 `````` ! assumes that the zeros at the beginning are passed through m1 res = GPL_zero_zi(m(1),y) `````` Luca committed Apr 09, 2019 483 484 `````` return end if `````` Luca committed Apr 09, 2019 485 `````` `````` 486 487 488 `````` ! has trailing zeroes? if(abs(z(k)) < zero ) then ! we remove them in flat form `````` Luca committed Apr 25, 2019 489 `````` z_flat = get_flattened_z(m,z) `````` 490 `````` res = G_flat(z_flat,y) `````` Luca committed May 22, 2019 491 `````` return `````` Luca committed Apr 25, 2019 492 493 `````` end if `````` Luca committed Apr 23, 2019 494 `````` ! need make convergent? `````` Luca Naterop committed Jul 08, 2019 495 `````` if(.not. all(abs(y) <= abs(z))) then `````` Luca committed May 22, 2019 496 497 `````` z_flat = get_flattened_z(m,z) res = G_flat(z_flat,y) `````` Luca committed Apr 25, 2019 498 `````` return `````` Luca committed Apr 09, 2019 499 `````` end if `````` ulrich_y committed Aug 28, 2019 500 `````` ! This is okay because in the MPLs the ieps doesn't matter anymore `````` ulrich_y committed Jul 09, 2019 501 `````` x(1) = y%c/z(1)%c `````` ulrich_y committed Jul 08, 2019 502 `````` do i = 2,k `````` ulrich_y committed Jul 09, 2019 503 `````` x(i) = z(i-1)%c/z(i)%c `````` Luca Naterop committed Apr 04, 2019 504 `````` end do `````` Luca committed May 20, 2019 505 506 507 `````` ! print*, 'computed using Li with ' ! print*, 'm = ', m ! print*, 'x = ', x `````` Luca committed Apr 25, 2019 508 `````` res = (-1)**k * MPL(m,x) `````` 509 `````` END FUNCTION G_condensed `````` Luca committed Apr 04, 2019 510 `````` `````` ulrich_y committed Jul 08, 2019 511 512 513 514 `````` FUNCTION G_SUPERFLATN(c0,n) integer, intent(in) :: n `````` ulrich_y committed Jul 09, 2019 515 `````` type(inum), intent(in) :: c0(n) `````` ulrich_y committed Jul 08, 2019 516 `````` complex(kind=prec) g_superflatn `````` ulrich_y committed Jul 09, 2019 517 `````` G_superflatn = G_flat(c0(1:n-1), c0(n)) `````` ulrich_y committed Jul 08, 2019 518 519 `````` END FUNCTION `````` ulrich_y committed Jul 19, 2019 520 521 522 `````` FUNCTION G_FLATr(Z_FLAT,Y) real(kind=prec), intent(in) :: z_flat(:), y complex(kind=prec) :: g_flatr `````` ulrich_y committed Aug 22, 2019 523 `````` g_flatr = G_flat(toinum(z_flat), toinum(y)) `````` ulrich_y committed Jul 19, 2019 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 committed Sep 27, 2019 529 `````` g_flatc = G_flat(toinum(z_flat), toinum(y)) `````` ulrich_y committed Jul 19, 2019 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 committed Apr 04, 2019 568 ``````END MODULE gpl_module `````` Luca Naterop committed Jul 07, 2019 569