fix transform to condensed notation, correct implementation of removal of trailing zeroes, and more

gpl_module.f90

@@ -6,7 +6,6 @@ MODULE gpl_module
 
 CONTAINS 
 
-
   RECURSIVE FUNCTION factorial(n) result(res)
     integer, intent(in) :: n
     integer :: res
@@ -37,7 +36,6 @@ CONTAINS
         GPL_has_convergent_series = .true.
       end if
     end if
-
   END FUNCTION GPL_has_convergent_series
 
   FUNCTION GPL_zero_zi(l,y)
@@ -46,15 +44,39 @@ CONTAINS
     complex(kind=prec) :: y, GPL_zero_zi
     print*, 'computed value using zi = 0'
     GPL_zero_zi = 1.0d0/factorial(l) * log(y) ** l
-
   END FUNCTION GPL_zero_zi
 
-  FUNCTION  G_with_flat_args(z_flat,y) result(res)
+  RECURSIVE FUNCTION G_flat(z_flat,y) result(res)
     ! Calls G function with flat arguments, that is, zeroes not passed through the m's. 
     complex(kind=prec) :: z_flat(:), y, res
-    complex(kind=prec), allocatable :: z(:)
-    integer :: m_prime(size(z_flat)), condensed_size
+    complex(kind=prec), allocatable :: z(:), s(:,:)
+    integer :: m_prime(size(z_flat)), condensed_size, kminusj, j, k, i
     integer, allocatable :: m(:)
+
+    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
     m_prime = get_condensed_m(z_flat)
     if(find_first_zero(m_prime) == -1) then
       condensed_size = size(m_prime)
@@ -65,40 +87,35 @@ CONTAINS
     allocate(z(condensed_size))
     m = m_prime(1:condensed_size)
     z = get_condensed_z(m,z_flat)
-    res = GPL(m,z,y,size(m))
+
+    res = G_condensed(m,z,y,size(m))
     deallocate(m)
     deallocate(z)
-  END  FUNCTION G_with_flat_args
+  END  FUNCTION G_flat
 
-  RECURSIVE FUNCTION GPL(m,z,y,k) result(res)
+  RECURSIVE FUNCTION G_condensed(m,z,y,k) result(res)
     ! computes the generalized polylogarithm G_{m1,..mk} (z1,...zk; y)
     ! assumes zero arguments expressed through the m's
     
-    integer :: m(:), k, i, kminusj
+    integer :: m(:), k, i
     complex(kind=prec) :: z(:), x(k), y, res, c(sum(m)+1,sum(m)+1), z_flat(sum(m)), a(sum(m)-1)
 
-    ! print*, 'z = ', abs(get_flattened_z(m,z))
+    print*, 'called G_condensed with args'
+    print*, 'm = ', m
+    print*, 'z = ', abs(z)
+
     ! are all z_i = 0 ? 
-    if(k == 1 .and. abs(z(1)) == 0) then
+    if(k == 1 .and. abs(z(1)) < zero) then
       ! assumes that the zeros at the beginning are passed through m1
       res = GPL_zero_zi(m(1),y)
       return
     end if
 
-    !  need to remove trailing  zeros?
-    if(abs(z(k)) == 0 ) then
-      print*, 'need to remove trailing zeros'   ! which we do in flat form
-      ! flatten z
+    ! has trailing zeroes?
+    if(abs(z(k)) < zero ) then
+      ! we remove them in flat form
       z_flat = get_flattened_z(m,z)
-      res = G_with_flat_args(z_flat,y)
-
-      ! a = z_flat(1: (size(z_flat)-1))
-      ! c = shuffle_with_zero(a)
-      ! res = G_with_flat_args(a,y)*log(y)
-      ! do  i = 2,k
-      !   res = res - G_with_flat_args(c(i,:),y)
-      ! end do
-      ! return
+      res = G_flat(z_flat,y)
     end  if
 
     ! need make convergent?
@@ -118,8 +135,7 @@ CONTAINS
     end do
     print*, 'computed using MPL'
     res = (-1)**k * MPL(m,x)
-
-  END FUNCTION GPL
+  END FUNCTION G_condensed
 
 END MODULE gpl_module
 

test.f90

@@ -67,7 +67,7 @@ CONTAINS
     character(len=*) :: test_id
 
     print*, '  ', 'testing GPL ', test_id, ' ...'
-    res = GPL(m,z,y,k)
+    res = G_condensed(m,z,y,k)
     call check(res,ref)
   end subroutine test_one_GPL
 
@@ -76,17 +76,18 @@ CONTAINS
     complex(kind=prec), parameter :: epsilon = 1E-14
     print*, 'doing GPL tests...'
     
-    ref = dcmplx(0.0819393734128676)
-    call test_one_GPL((/ 1,1 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.1')
+    ! ref = dcmplx(0.0819393734128676)
+    ! call test_one_GPL((/ 1,1 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.1')
     
-    ref = dcmplx(0.01592795952537145)
-    call test_one_GPL((/ 3,2 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.2')
+    ! ref = dcmplx(0.01592795952537145)
+    ! call test_one_GPL((/ 3,2 /),cmplx((/ 1.3d0, 1.1d0 /)),cmplx(0.4),2,ref,'2.2')
     
-    ref = dcmplx(0.0020332632172573974)
-    call test_one_GPL((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3')
+    ! ref = dcmplx(0.0020332632172573974)
+    ! call test_one_GPL((/ 4 /),cmplx((/ 0 /)),cmplx(1.6),1,ref,'2.3')
     
     ref = dcmplx(0.0020332632172573974)
-    call test_one_GPL((/ 1,1 /),cmplx((/ 1.7,0.0 /)),cmplx(1.1),2,ref,'2.4')
+    ref = G_flat(cmplx((/1.7,0.0,0.0/)),cmplx(1.1))
+    ! call test_one_GPL((/1,1,1/),cmplx((/ 0.0,1.7,0.0 /)),cmplx(1.1),3,ref,'2.4')
 
   end subroutine do_GPL_tests
 

utils.f90

@@ -12,7 +12,7 @@
 MODULE utils
   implicit none
   integer, parameter :: prec = selected_real_kind(15,32)  
-
+  real :: zero = 1e-15
   ! logical :: print_enabled = .true.
   ! logical :: warnings_enabled = .true.
 
@@ -20,11 +20,12 @@ CONTAINS
 
   FUNCTION  get_condensed_m(z) result(m)
     ! returns condensed m where the ones not needed are filled with 0
-    complex(kind=prec) :: z(:), m(size(z))
-    integer :: pos = 1, i 
+    complex(kind=prec), intent(in) :: z(:)
+    integer :: m(size(z)), pos, i 
     m = 1
+    pos = 1
     do i = 1,size(z)
-      if(z(i) == 0) then
+      if(abs(z(i)) < zero) then
         if(i == size(z)) then
           pos = pos + 1
         else 
@@ -65,7 +66,7 @@ CONTAINS
     integer :: res, i
     res = 0
     do i = size(z), 1, -1
-      if( z(i) == 0 ) then
+      if( abs(z(i)) < zero ) then
         res = res + 1
       else
         exit
@@ -91,10 +92,16 @@ CONTAINS
     integer :: s(2), i
     s = shape(m)
     do i = 1,s(1)
-      print*, m(i,:)
+      print*, abs(m(i,:))
     end do
   END SUBROUTINE print_as_matrix
 
+  FUNCTION zero_array(n) result(res)
+    integer :: n
+    complex(kind=prec) :: res(n)
+    res = 0
+  END FUNCTION zero_array
+
   FUNCTION shuffle_with_zero(a) result(res)
     ! rows of result are shuffles of a with 0
     complex :: a(:)
@@ -130,35 +137,21 @@ END MODULE utils
 ! PROGRAM test
 !   use  utils
 !   implicit none
-
-
-!   ! complex(kind=prec), dimension(5) :: a = cmplx((/1,2,3/))
-!   ! complex(kind=prec) :: z_flat(2)
-!   ! complex(kind=prec), allocatable :: z(:)
-!   ! integer :: m_prime(2), condensed_size
-!   ! integer, allocatable :: m(:)
-!   ! complex(kind=prec) :: b(size(a)+1,size(a)+1)
-
-!   ! ! ! test shuffling
-!   ! ! b = 1
-!   ! ! b = shuffle_with_zero(a)
-!   ! ! call print_as_matrix(b)
-
-!   ! ! test condensing
-!   ! z_flat = cmplx((/4,0/))
-!   ! 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)
-!   ! z_flat = get_flattened_z(m,z)
-!   ! deallocate(m)
-!   ! deallocate(z)
+  
+!   complex(kind=prec) :: z(4)
+!   integer :: m_prime(4), condensed_size
+!   z = cmplx((/0.0,1.7,0.0,0.0/))
+
+!   ! transform to condensed notation
+!   m_prime = get_condensed_m(z)
+!   print*, abs(z)
+!   m_prime = get_condensed_m(z)
+!   print*, abs(z)
+!   if(find_first_zero(m_prime) == -1) then
+!     condensed_size = size(m_prime)
+!   else
+!     condensed_size = find_first_zero(m_prime)-1 
+!   end if
+!   print*, condensed_size
 
 ! END  PROGRAM test
-