Added overflow control to naive_polylog

src/gpl_module.f90

@@ -67,7 +67,7 @@ CONTAINS
     if(.not. present(y)) print*, 'G(', abs(z_flat), ')'
   END SUBROUTINE print_G
 
-  FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) result(res)
+  RECURSIVE FUNCTION remove_sr_from_last_place_in_PI(a,y2,m,p) result(res)
     ! here what is passed is not the full a vector, only a1, ..., ak without the trailing zeroes
     integer :: m, i, j, n
     complex(kind=prec) :: a(:), y2, s(m), p(:), res
@@ -127,7 +127,7 @@ CONTAINS
       return
     end if
   
-    a = g(1:size(p)-1)
+    a = g(1:size(g)-1)
     y2 = g(size(g)) 
     m = size(g)  ! actually, size(g)-1+1
 

src/maths_functions.f90

@@ -11,10 +11,13 @@ CONTAINS
     integer :: m
     complex(kind=prec) :: x, res
     integer :: i,n
-    integer, allocatable :: j(:)
     n = 1000
-    j = (/(i, i=1,n,1)/) 
-    res = sum(x**j / j**m)
+    res=0.
+    do i=1,n
+      if(i**m.lt.0) return ! roll over
+      if(abs(x**i).lt.1.e-250) return
+      res = res+x**i/i**m
+    enddo
   END FUNCTION naive_polylog
 
   FUNCTION Li2(x)