gpl_module.f90 17.4 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 Jul 09, 2019 55 `````` RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) 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 Jul 08, 2019 58 59 `````` type(inum) :: a(:), y2, s(m), p(:) complex(kind=prec) :: res `````` ulrich_y committed Jul 08, 2019 60 `````` type(inum) :: alpha(product((/(i,i=1,size(a)+size(s))/))/ & `````` Luca Naterop committed Jul 05, 2019 61 62 63 `````` (product((/(i,i=1,size(a))/))*product((/(i,i=1,size(s))/))), & size(a) + size(s)) `````` ulrich_y committed Jul 08, 2019 64 `````` s = [zeroes(m-1),marker] `````` Luca Naterop committed Jul 05, 2019 65 `````` alpha = shuffle_product(a,s) `````` ulrich_y committed Jul 10, 2019 66 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 67 68 69 `````` if(verb >= 50) then print*, 'mapping to ' call print_G(a,y2) `````` Luca Naterop committed Jul 08, 2019 70 `````` print*, 'PI with p=',real(p),'i=',m,'g =',real([zeroes(m-1),y2]) `````` Luca Naterop committed Jul 05, 2019 71 `````` end if `````` ulrich_y committed Jul 10, 2019 72 ``````#endif `````` ulrich_y committed Jul 10, 2019 73 `````` res = G_flat(a,y2)*pending_integral(p,m,[zeroes(m-1),y2]) `````` ulrich_y committed Jul 10, 2019 74 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 75 `````` if(verb >= 50) print*, 'also mapping to' `````` ulrich_y committed Jul 10, 2019 76 ``````#endif `````` Luca Naterop committed Jul 05, 2019 77 78 `````` do j = 2,size(alpha, 1) ! find location of s_r `````` ulrich_y committed Jul 08, 2019 79 `````` n = find_marker(alpha(j,:)) `````` ulrich_y committed Jul 10, 2019 80 ``````#ifdef DEBUG `````` ulrich_y committed Jul 08, 2019 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 committed Jul 10, 2019 83 ``````#endif `````` Luca Naterop committed Jul 05, 2019 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 `````` Luca Naterop committed May 16, 2019 88 `````` RECURSIVE FUNCTION pending_integral(p,i,g) result(res) `````` Luca committed May 21, 2019 89 `````` ! evaluates a pending integral by reducing it to simpler ones and g functions `````` ulrich_y committed Jul 08, 2019 90 91 92 `````` 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 93 `````` integer :: i, m `````` Luca committed May 21, 2019 94 `````` res = 0 `````` Luca committed May 23, 2019 95 `````` `````` ulrich_y committed Jul 10, 2019 96 ``````#ifdef DEBUG `````` ulrich_y committed Jul 08, 2019 97 `````` if(verb >= 30) print*, 'evaluating PI with p=',real(p),'i=',real(i),'g =',real(g) `````` ulrich_y committed Jul 10, 2019 98 ``````#endif `````` Luca Naterop committed Jul 05, 2019 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 committed Jul 10, 2019 105 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 106 `````` if(verb >= 30) print*, 'only integrals in front, we get G-function' `````` ulrich_y committed Jul 10, 2019 107 ``````#endif `````` Luca Naterop committed Jul 05, 2019 108 109 110 `````` res = G_flat(b,y1) return end if `````` Luca committed Jun 04, 2019 111 `````` `````` Luca committed May 23, 2019 112 `````` ! if integration variable at end -> we gat a G function `````` Luca committed May 21, 2019 113 `````` if(i == size(g)+1) then `````` ulrich_y committed Jul 10, 2019 114 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 115 `````` if(verb >= 30) print*, 'is just a G-function' `````` ulrich_y committed Jul 10, 2019 116 ``````#endif `````` Luca committed May 21, 2019 117 118 119 `````` res = G_flat([p(2:size(p)),g], p(1)) return end if `````` Luca committed Jun 04, 2019 120 `````` `````` Luca Naterop committed Jul 05, 2019 121 `````` `````` Luca Naterop committed Aug 28, 2019 122 `````` ! if depth one and m = 1 use (23) `````` Luca committed May 21, 2019 123 `````` if(size(g) == 1) then `````` ulrich_y committed Jul 10, 2019 124 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 125 `````` if(verb >= 30) print*, 'case depth one and m = 1' `````` ulrich_y committed Jul 10, 2019 126 ``````#endif `````` ulrich_y committed Jul 08, 2019 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 committed Aug 28, 2019 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_) `````` Luca committed May 21, 2019 131 132 `````` return end if `````` Luca committed Jun 04, 2019 133 `````` `````` ulrich_y committed Jul 09, 2019 134 `````` a = g(1:size(g)-1) `````` Luca Naterop committed Jul 05, 2019 135 136 137 `````` y2 = g(size(g)) m = size(g) ! actually, size(g)-1+1 `````` Luca Naterop committed Aug 28, 2019 138 `````` ! if depth one and m > 1 use (25) `````` Luca Naterop committed Jul 05, 2019 139 `````` if(all( abs( g(1:size(g)-1) ) < zero)) then `````` ulrich_y committed Jul 10, 2019 140 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 141 142 143 144 `````` 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 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 committed Jul 08, 2019 148 `````` print*, 'PI with p=',tocmplx([p,izero]),'i=',m-1,'g =',tocmplx(zeroes(0)) `````` Luca Naterop committed Jul 05, 2019 149 `````` end if `````` ulrich_y committed Jul 10, 2019 150 ``````#endif `````` Luca Naterop committed Jul 05, 2019 151 `````` res = -zeta(m)*pending_integral(p,0,zeroes(0)) & `````` ulrich_y committed Jul 08, 2019 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 committed Jun 04, 2019 154 155 156 `````` return end if `````` Luca Naterop committed Aug 28, 2019 157 `````` ! case of higher depth, s_r at beginning, use (21) `````` Luca committed May 23, 2019 158 `````` if(i == 1) then `````` ulrich_y committed Jul 10, 2019 159 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 05, 2019 160 `````` if(verb >= 30) print*, 'case higher depth, sr at beginning' `````` Luca committed Jun 04, 2019 161 `````` `````` Luca Naterop committed Jul 05, 2019 162 `````` if(verb >= 50) then `````` Luca Naterop committed Aug 28, 2019 163 `````` print*, 'map to (using 21)' `````` ulrich_y committed Jul 08, 2019 164 `````` print*, 'PI with p=',real(p),'i=',0,'g =',tocmplx(zeroes(0) ) `````` ulrich_y committed Jul 08, 2019 165 `````` call print_G([izero,a],y2) `````` ulrich_y committed Jul 08, 2019 166 `````` print*, 'PI with p=',real([p,y2]),'i=',0,'g =',tocmplx(zeroes(0) ) `````` Luca Naterop committed Jul 05, 2019 167 `````` call print_G(a,y2) `````` Luca Naterop committed Jul 08, 2019 168 `````` print*, 'PI with p=',real([p,a(1)]),'i=',1,'g =',real([a(2:size(a)),y2]) `````` ulrich_y committed Jul 08, 2019 169 `````` print*, 'PI with p=',real([p,a(1)]),'i=',0,'g =',tocmplx(zeroes(0)) `````` Luca Naterop committed Jul 05, 2019 170 171 `````` call print_G(a,y2) end if `````` ulrich_y committed Jul 10, 2019 172 ``````#endif `````` Luca Naterop committed Jul 05, 2019 173 `````` `````` ulrich_y committed Jul 08, 2019 174 `````` res = pending_integral(p,0,zeroes(0)) * G_flat([izero,a],y2) & `````` Luca Naterop committed Jul 05, 2019 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) `````` Luca committed May 23, 2019 178 179 `````` return end if `````` Luca committed Jun 04, 2019 180 `````` `````` Luca Naterop committed Aug 28, 2019 181 `````` ! case higher depth, s_r at the end, use shuffle algebra to remove from last place `````` Luca committed Jun 04, 2019 182 `````` if(i == size(g)) then `````` ulrich_y committed Jul 10, 2019 183 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 184 `````` if(verb >= 30) print*, 's_r at the end under PI, need to shuffle' `````` ulrich_y committed Jul 10, 2019 185 ``````#endif `````` Luca Naterop committed Jul 05, 2019 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 committed Jun 04, 2019 188 189 190 `````` return end if `````` Luca Naterop committed Aug 28, 2019 191 `````` ! case higher depth, s_r in middle, use (22) `````` ulrich_y committed Jul 10, 2019 192 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 193 `````` if(verb >= 30) print*, 's_r in the middle under PI' `````` ulrich_y committed Jul 10, 2019 194 ``````#endif `````` Luca Naterop committed Jul 07, 2019 195 `````` `````` Luca Naterop committed Aug 28, 2019 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 committed May 15, 2019 201 202 `````` END FUNCTION pending_integral `````` Luca Naterop committed Jul 08, 2019 203 `````` RECURSIVE FUNCTION remove_sr_from_last_place_in_G(a,y2,m,sr) result(res) `````` ulrich_y committed Jul 08, 2019 204 205 `````` type(inum) :: a(:), sr, y2 complex(kind=prec) :: res `````` Luca committed May 21, 2019 206 `````` integer :: m,i,j `````` ulrich_y committed Jul 08, 2019 207 `````` type(inum) :: alpha(product((/(i,i=1,size(a)+m)/))/ & `````` Luca committed May 21, 2019 208 209 `````` (product((/(i,i=1,size(a))/))*product((/(i,i=1,m)/))), & size(a) + m) `````` Luca committed Jun 04, 2019 210 211 `````` 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 212 213 214 `````` do j = 2,size(alpha,1) res = res - G_flat(alpha(j,:),y2) end do `````` Luca Naterop committed Jul 05, 2019 215 `````` END FUNCTION remove_sr_from_last_place_in_G `````` Luca committed May 21, 2019 216 `````` `````` Luca Naterop committed Jul 08, 2019 217 `````` RECURSIVE FUNCTION make_convergent(a,y2) result(res) `````` Luca Naterop committed Aug 28, 2019 218 `````` ! goes from G-functions to pending integrals and simpler expressions `````` Luca committed Jun 04, 2019 219 `````` `````` ulrich_y committed Jul 08, 2019 220 221 `````` type(inum) :: a(:), y2, sr complex(kind=prec) :: res `````` Luca committed May 21, 2019 222 `````` integer :: i, mminus1 `````` Luca committed May 20, 2019 223 `````` `````` Luca committed May 15, 2019 224 `````` res = 0 `````` Luca Naterop committed May 16, 2019 225 226 `````` i = min_index(abs(a)) sr = a(i) `````` Luca committed Jun 04, 2019 227 `````` `````` Luca Naterop committed Jul 07, 2019 228 `````` if(i == size(a)) then `````` Luca Naterop committed Aug 28, 2019 229 `````` ! sr at the end, thus shuffle `````` ulrich_y committed Jul 10, 2019 230 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 231 `````` if(verb >= 30) print*, 'sr at the end' `````` ulrich_y committed Jul 10, 2019 232 ``````#endif `````` Luca Naterop committed Jul 07, 2019 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 committed Jun 04, 2019 237 `````` `````` Luca Naterop committed May 16, 2019 238 `````` if(i == 1) then `````` Luca Naterop committed Aug 28, 2019 239 `````` !s_r at beginning, thus use (18) `````` Luca committed Jun 04, 2019 240 `````` `````` ulrich_y committed Jul 10, 2019 241 ``````#ifdef DEBUG `````` Luca committed Jun 04, 2019 242 243 244 `````` if(verb >= 30) then print*, '--------------------------------------------------' print*, 'sr at beginning, map to: ' `````` ulrich_y committed Jul 08, 2019 245 `````` call print_G([izero, a(i+1:size(a))], y2) `````` Luca committed Jun 04, 2019 246 247 `````` call print_G([y2], sr) call print_G(a(i+1:size(a)), y2) `````` Luca Naterop committed Jul 08, 2019 248 `````` print*, 'PI with p=',real([sr, a(i+1)]),'i=',i,'g =', real([a(i+2:size(a)), y2]) `````` Luca committed Jun 04, 2019 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 committed Jul 10, 2019 253 ``````#endif `````` Luca committed Jun 04, 2019 254 `````` `````` ulrich_y committed Jul 08, 2019 255 `````` res = G_flat([izero, a(i+1:size(a))], y2) & `````` Luca Naterop committed May 16, 2019 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 committed Jul 08, 2019 258 `````` - G_flat([a(i+1)],sr) * G_flat(a(i+1:size(a)), y2) `````` Luca committed May 15, 2019 259 260 `````` return end if `````` Luca committed May 21, 2019 261 `````` `````` Luca Naterop committed Aug 28, 2019 262 `````` ! so s_r in middle, use (19) `````` ulrich_y committed Jul 10, 2019 263 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 264 265 266 `````` if(verb >= 30) then print*, '--------------------------------------------------' print*, 'sr in the middle, map to: ' `````` ulrich_y committed Jul 08, 2019 267 `````` call print_G([a(1:i-1),izero,a(i+1:size(a))],y2) `````` Luca Naterop committed Jul 08, 2019 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]) `````` Luca committed May 23, 2019 269 270 `````` 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 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]) `````` Luca committed May 23, 2019 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 committed Jul 10, 2019 276 ``````#endif `````` Luca committed May 23, 2019 277 `````` `````` ulrich_y committed Jul 08, 2019 278 `````` res = G_flat([a(1:i-1),izero,a(i+1:size(a))],y2) & `````` Luca committed May 21, 2019 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 committed Jun 04, 2019 283 `````` END FUNCTION make_convergent `````` Luca committed May 15, 2019 284 `````` `````` ulrich_y committed Jul 08, 2019 285 `````` RECURSIVE FUNCTION improve_convergence(z) result(res) `````` Luca Naterop committed Jul 07, 2019 286 `````` ! improves the convergence by applying the Hoelder convolution to G(z1,...zk,1) `````` ulrich_y committed Jul 08, 2019 287 288 `````` type(inum) :: z(:),oneminusz(size(z)) complex(kind=prec) :: res `````` Luca Naterop committed Jul 07, 2019 289 290 `````` complex(kind=prec), parameter :: p = 2.0 integer :: k, j `````` ulrich_y committed Jul 10, 2019 291 ``````#ifdef DEBUG `````` Luca Naterop committed Jul 07, 2019 292 `````` if(verb >= 30) print*, 'requires Hoelder convolution' `````` ulrich_y committed Jul 10, 2019 293 ``````#endif `````` ulrich_y committed Jul 09, 2019 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 committed Jul 09, 2019 298 `````` do j=1,size(z) `````` ulrich_y committed Jul 09, 2019 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 committed Jul 09, 2019 304 `````` enddo `````` Luca Naterop committed Jul 07, 2019 305 306 `````` k = size(z) `````` ulrich_y committed Jul 09, 2019 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 committed Jul 07, 2019 309 `````` do j = 1,k-1 `````` ulrich_y committed Jul 09, 2019 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 committed Jul 07, 2019 311 312 313 `````` end do END FUNCTION improve_convergence `````` 314 `````` RECURSIVE FUNCTION G_flat(z_flat,y) result(res) `````` Luca Naterop committed Apr 29, 2019 315 `````` ! Calls G function with flat arguments, that is, zeroes not passed through the m's. `````` ulrich_y committed Jul 09, 2019 316 `````` type(inum) :: z_flat(:), y, znorm(size(z_flat)) `````` ulrich_y committed Jul 08, 2019 317 318 `````` complex(kind=prec) :: res type(inum), allocatable :: s(:,:), z(:) `````` Luca committed May 20, 2019 319 `````` integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i, m_1 `````` Luca committed Apr 25, 2019 320 `````` integer, allocatable :: m(:) `````` Luca committed May 20, 2019 321 `````` logical :: is_depth_one `````` 322 `````` `````` ulrich_y committed Jul 10, 2019 323 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 324 `````` if(verb >= 50) call print_G(z_flat,y) `````` ulrich_y committed Jul 10, 2019 325 ``````#endif `````` Luca committed May 15, 2019 326 `````` `````` Luca Naterop committed Jul 08, 2019 327 `````` `````` Luca Naterop committed Jul 08, 2019 328 `````` if(size(z_flat) == 1) then `````` ulrich_y committed Jul 09, 2019 329 `````` if( abs(z_flat(1)%c - y%c) <= zero ) then `````` Luca Naterop committed Jul 08, 2019 330 331 332 333 334 `````` res = 0 return end if end if `````` Luca Naterop committed Jul 05, 2019 335 `````` ! add small imaginary part if not there `````` Luca Naterop committed Jul 05, 2019 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 committed May 21, 2019 340 `````` `````` Luca committed May 21, 2019 341 `````` ! is just a logarithm? `````` Luca committed May 21, 2019 342 `````` if(all(abs(z_flat) < zero)) then `````` ulrich_y committed Jul 10, 2019 343 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 344 `````` if(verb >= 70) print*, 'all z are zero' `````` ulrich_y committed Jul 10, 2019 345 ``````#endif `````` ulrich_y committed Jul 09, 2019 346 `````` res = gpl_zero_zi(size(z_flat),y) `````` Luca committed May 21, 2019 347 348 `````` return end if `````` Luca committed May 21, 2019 349 `````` if(size(z_flat) == 1) then `````` ulrich_y committed Jul 10, 2019 350 ``````#ifdef DEBUG `````` Luca committed May 23, 2019 351 `````` if(verb >= 70) print*, 'is just a logarithm' `````` ulrich_y committed Jul 10, 2019 352 ``````#endif `````` ulrich_y committed Jul 09, 2019 353 `````` res = plog1(y,z_flat(1)) ! log((z_flat(1) - y)/z_flat(1)) `````` Luca committed May 21, 2019 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 committed Jul 10, 2019 363 ``````#ifdef DEBUG `````` Luca committed May 22, 2019 364 `````` if(verb >= 70) print*, 'is just a polylog' `````` ulrich_y committed Jul 10, 2019 365 ``````#endif `````` ulrich_y committed Jul 09, 2019 366 `````` res = -polylog(m_1, y, z_flat(m_1))!-polylog(m_1,y/z_flat(m_1)) `````` Luca committed May 21, 2019 367 368 369 `````` return end if `````` Luca committed May 07, 2019 370 `````` ! need remove trailing zeroes? `````` 371 372 373 `````` k = size(z_flat) kminusj = find_amount_trailing_zeros(z_flat) j = k - kminusj `````` ulrich_y committed Aug 19, 2019 374 `````` if(kminusj > 0) then `````` ulrich_y committed Jul 10, 2019 375 ``````#ifdef DEBUG `````` Luca committed May 22, 2019 376 `````` if(verb >= 30) print*, 'need to remove trailing zeroes' `````` ulrich_y committed Jul 10, 2019 377 ``````#endif `````` 378 379 `````` allocate(s(j,j)) s = shuffle_with_zero(z_flat(1:j-1)) `````` ulrich_y committed Jul 09, 2019 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 committed Jun 04, 2019 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 committed May 21, 2019 389 390 `````` ! need make convergent? if(.not. is_convergent(z_flat,y)) then `````` ulrich_y committed Jul 10, 2019 391 ``````#ifdef DEBUG `````` Luca committed May 22, 2019 392 `````` if(verb >= 10) print*, 'need to make convergent' `````` ulrich_y committed Jul 10, 2019 393 ``````#endif `````` Luca committed Jun 04, 2019 394 `````` res = make_convergent(z_flat, y) `````` Luca committed May 21, 2019 395 396 397 `````` return end if `````` Luca Naterop committed Jul 07, 2019 398 `````` ! requires Hoelder convolution? `````` ulrich_y committed Jul 10, 2019 399 `````` 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 400 `````` ! Here we *assume* that y is positive and doesn't mess up the `````` ulrich_y committed Jul 09, 2019 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 committed Jul 07, 2019 406 407 408 `````` return end if `````` Luca committed May 15, 2019 409 `````` ! thus it is convergent, and has no trailing zeros `````` Luca Naterop committed Jul 10, 2019 410 `````` ! -> evaluate in condensed notation which will give series representation `````` Luca committed Apr 25, 2019 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)) `````` Luca committed Apr 25, 2019 422 423 `````` deallocate(m) deallocate(z) `````` Luca Naterop committed May 16, 2019 424 `````` END FUNCTION G_flat `````` Luca committed Apr 25, 2019 425 `````` `````` Luca committed May 22, 2019 426 427 428 `````` FUNCTION G_superflat(g) result(res) ! simpler notation for flat evaluation complex(kind=prec) :: g(:), res `````` ulrich_y committed Jul 08, 2019 429 `````` res = G_flat(toinum(g(1:size(g)-1)), inum(g(size(g)),di0)) `````` Luca committed May 22, 2019 430 431 `````` END FUNCTION G_superflat `````` Luca committed May 23, 2019 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 committed Aug 22, 2019 436 `````` res = G_flat(toinum(g(1:size(g)-1)), toinum(g(size(g)))) `````` Luca committed May 23, 2019 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 committed Aug 22, 2019 443 `````` res = G_flat(toinum(real(g(1:size(g)-1),kind=prec)), toinum(real(g(size(g)),kind=prec))) `````` Luca committed May 23, 2019 444 445 `````` END FUNCTION G_int `````` 446 `````` RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res) `````` Luca Naterop committed Apr 04, 2019 447 `````` ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) `````` Luca committed Apr 09, 2019 448 `````` ! assumes zero arguments expressed through the m's `````` Luca committed Apr 23, 2019 449 `````` `````` 450 `````` integer :: m(:), k, i `````` ulrich_y committed Jul 09, 2019 451 452 `````` type(inum) :: z(:), y, z_flat(sum(m)) complex(kind=prec) :: res, x(k) `````` Luca committed Apr 25, 2019 453 `````` `````` Luca committed May 08, 2019 454 455 456 `````` ! print*, 'called G_condensed with args' ! print*, 'm = ', m ! print*, 'z = ', abs(z) `````` 457 `````` `````` Luca committed Apr 09, 2019 458 `````` ! are all z_i = 0 ? `````` 459 `````` if(k == 1 .and. abs(z(1)) < zero) then `````` Luca committed Apr 25, 2019 460 461 `````` ! assumes that the zeros at the beginning are passed through m1 res = GPL_zero_zi(m(1),y) `````` Luca committed Apr 09, 2019 462 463 `````` return end if `````` Luca committed Apr 09, 2019 464 `````` `````` 465 466 467 `````` ! has trailing zeroes? if(abs(z(k)) < zero ) then ! we remove them in flat form `````` Luca committed Apr 25, 2019 468 `````` z_flat = get_flattened_z(m,z) `````` 469 `````` res = G_flat(z_flat,y) `````` Luca committed May 22, 2019 470 `````` return `````` Luca committed Apr 25, 2019 471 472 `````` end if `````` Luca committed Apr 23, 2019 473 `````` ! need make convergent? `````` Luca Naterop committed Jul 08, 2019 474 `````` if(.not. all(abs(y) <= abs(z))) then `````` Luca committed May 22, 2019 475 476 `````` z_flat = get_flattened_z(m,z) res = G_flat(z_flat,y) `````` Luca committed Apr 25, 2019 477 `````` return `````` Luca committed Apr 09, 2019 478 `````` end if `````` ulrich_y committed Aug 28, 2019 479 `````` ! This is okay because in the MPLs the ieps doesn't matter anymore `````` ulrich_y committed Jul 09, 2019 480 `````` x(1) = y%c/z(1)%c `````` ulrich_y committed Jul 08, 2019 481 `````` do i = 2,k `````` ulrich_y committed Jul 09, 2019 482 `````` x(i) = z(i-1)%c/z(i)%c `````` Luca Naterop committed Apr 04, 2019 483 `````` end do `````` Luca committed May 20, 2019 484 485 486 `````` ! print*, 'computed using Li with ' ! print*, 'm = ', m ! print*, 'x = ', x `````` Luca committed Apr 25, 2019 487 `````` res = (-1)**k * MPL(m,x) `````` 488 `````` END FUNCTION G_condensed `````` Luca committed Apr 04, 2019 489 `````` `````` ulrich_y committed Jul 08, 2019 490 491 492 493 `````` FUNCTION G_SUPERFLATN(c0,n) integer, intent(in) :: n `````` ulrich_y committed Jul 09, 2019 494 `````` type(inum), intent(in) :: c0(n) `````` ulrich_y committed Jul 08, 2019 495 `````` complex(kind=prec) g_superflatn `````` ulrich_y committed Jul 09, 2019 496 `````` G_superflatn = G_flat(c0(1:n-1), c0(n)) `````` ulrich_y committed Jul 08, 2019 497 498 `````` END FUNCTION `````` ulrich_y committed Jul 19, 2019 499 500 501 `````` 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 502 `````` g_flatr = G_flat(toinum(z_flat), toinum(y)) `````` ulrich_y committed Jul 19, 2019 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 committed Apr 04, 2019 547 ``````END MODULE gpl_module `````` Luca Naterop committed Jul 07, 2019 548