gpl_module.f90 4 KB
 Luca Naterop committed Apr 04, 2019 1 `````` `````` Luca committed Apr 04, 2019 2 3 ``````MODULE gpl_module use mpl_module `````` Luca committed Apr 25, 2019 4 `````` use utils `````` Luca committed Apr 04, 2019 5 6 7 8 `````` implicit none CONTAINS `````` Luca committed Apr 09, 2019 9 10 11 12 13 14 `````` RECURSIVE FUNCTION factorial(n) result(res) integer, intent(in) :: n integer :: res res = merge(1,n*factorial(n-1),n==1) END FUNCTION factorial `````` Luca committed Apr 09, 2019 15 16 17 18 19 20 21 22 23 `````` FUNCTION zeta(n) real(kind=prec) :: values(9), zeta integer :: n values = (/1.6449340668482262, 1.2020569031595942, 1.0823232337111381, & 1.03692775514337, 1.0173430619844488, 1.008349277381923, & 1.0040773561979441, 1.0020083928260821, 1.000994575127818/) zeta = values(n-1) END FUNCTION zeta `````` Luca committed Apr 09, 2019 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 `````` FUNCTION GPL_has_convergent_series(m,z,y,k) ! tests if GPL has a convergent series representation integer :: m(:), k complex(kind=prec) :: z(:), y logical :: GPL_has_convergent_series GPL_has_convergent_series = .false. if(all(abs(y) <= abs(z))) then if(m(1) == 1) then GPL_has_convergent_series = (y/z(1) /= 1) else GPL_has_convergent_series = .true. end if end if END FUNCTION GPL_has_convergent_series `````` Luca committed Apr 09, 2019 41 42 43 44 `````` FUNCTION GPL_zero_zi(l,y) ! used to compute the value of GPL when all zi are zero integer :: l complex(kind=prec) :: y, GPL_zero_zi `````` Luca Naterop committed Apr 29, 2019 45 `````` print*, 'computed value using zi = 0' `````` Luca committed Apr 09, 2019 46 47 48 `````` GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l END FUNCTION GPL_zero_zi `````` 49 `````` RECURSIVE FUNCTION G_flat(z_flat,y) result(res) `````` Luca Naterop committed Apr 29, 2019 50 `````` ! Calls G function with flat arguments, that is, zeroes not passed through the m's. `````` Luca committed Apr 25, 2019 51 `````` complex(kind=prec) :: z_flat(:), y, res `````` 52 53 `````` complex(kind=prec), allocatable :: z(:), s(:,:) integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i `````` Luca committed Apr 25, 2019 54 `````` integer, allocatable :: m(:) `````` 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 `````` print*, 'G_flat called with args', abs(z_flat) ! remove trailing zeroes k = size(z_flat) kminusj = find_amount_trailing_zeros(z_flat) j = k - kminusj if(all(abs(z_flat) < zero)) then print*, 'all are zero', abs(z_flat) res = GPL_zero_zi(k,y) return else if(kminusj > 0) then print*, 'we have',kminusj,'trailing zeroes' allocate(s(j,j)) s = shuffle_with_zero(z_flat(1:j-1)) res = log(y)*G_flat(z_flat(1:size(z_flat)-1),y) do i = 1,size(s,1) res = res - G_flat([s(i,:),z_flat(j),zero_array(kminusj-1)], y) end do res = res / kminusj deallocate(s) return end if ! transform to condensed notation `````` Luca committed Apr 25, 2019 80 81 82 83 84 85 86 87 88 89 `````` 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) `````` 90 91 `````` res = G_condensed(m,z,y,size(m)) `````` Luca committed Apr 25, 2019 92 93 `````` deallocate(m) deallocate(z) `````` 94 `````` END FUNCTION G_flat `````` Luca committed Apr 25, 2019 95 `````` `````` 96 `````` RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res) `````` Luca Naterop committed Apr 04, 2019 97 `````` ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y) `````` Luca committed Apr 09, 2019 98 `````` ! assumes zero arguments expressed through the m's `````` Luca committed Apr 23, 2019 99 `````` `````` 100 `````` integer :: m(:), k, i `````` Luca committed Apr 25, 2019 101 102 `````` complex(kind=prec) :: z(:), x(k), y, res, c(sum(m)+1,sum(m)+1), z_flat(sum(m)), a(sum(m)-1) `````` 103 104 105 106 `````` print*, 'called G_condensed with args' print*, 'm = ', m print*, 'z = ', abs(z) `````` Luca committed Apr 09, 2019 107 `````` ! are all z_i = 0 ? `````` 108 `````` if(k == 1 .and. abs(z(1)) < zero) then `````` Luca committed Apr 25, 2019 109 110 `````` ! assumes that the zeros at the beginning are passed through m1 res = GPL_zero_zi(m(1),y) `````` Luca committed Apr 09, 2019 111 112 `````` return end if `````` Luca committed Apr 09, 2019 113 `````` `````` 114 115 116 `````` ! has trailing zeroes? if(abs(z(k)) < zero ) then ! we remove them in flat form `````` Luca committed Apr 25, 2019 117 `````` z_flat = get_flattened_z(m,z) `````` 118 `````` res = G_flat(z_flat,y) `````` Luca committed Apr 25, 2019 119 120 `````` end if `````` Luca committed Apr 23, 2019 121 `````` ! need make convergent? `````` Luca committed Apr 09, 2019 122 `````` if(.not. GPL_has_convergent_series(m,z,y,k)) then `````` Luca committed Apr 23, 2019 123 124 125 126 127 128 `````` print*, 'need to make convergent' if(k == 1 .and. m(1) == 1) then print*, 'now we use the easy case. ha. ha.' else print*, ' ', 'does not have convergent series representation' end if `````` Luca committed Apr 25, 2019 129 130 `````` res = 0 return `````` Luca committed Apr 09, 2019 131 132 `````` end if `````` Luca Naterop committed Apr 04, 2019 133 134 135 `````` do i = 1,k x(i) = merge(y/z(1), z(i-1)/z(i),i == 1) end do `````` Luca committed Apr 25, 2019 136 137 `````` print*, 'computed using MPL' res = (-1)**k * MPL(m,x) `````` 138 `````` END FUNCTION G_condensed `````` Luca committed Apr 04, 2019 139 140 `````` END MODULE gpl_module `````` Luca Naterop committed Apr 04, 2019 141 142 `````` ``````