Commit 13945806 by ulrich_y

### Added ikin for integers prone to roll over

parent 991385d2
 ... @@ -4,6 +4,7 @@ MODULE globals ... @@ -4,6 +4,7 @@ MODULE globals integer, parameter :: prec = selected_real_kind(15,32) integer, parameter :: prec = selected_real_kind(15,32) integer, parameter :: ikin = 4 real(kind=prec), parameter :: zero = 10._prec**(-precision(1._prec)) ! values smaller than this count as zero real(kind=prec), parameter :: zero = 10._prec**(-precision(1._prec)) ! values smaller than this count as zero real(kind=prec), parameter :: pi = 3.1415926535897932384626433832795028841971693993751_prec real(kind=prec), parameter :: pi = 3.1415926535897932384626433832795028841971693993751_prec ... ...
 ... @@ -26,7 +26,7 @@ CONTAINS ... @@ -26,7 +26,7 @@ CONTAINS ! Computes the classical polylogarithm Li_m(x) using series representation up to order n ! Computes the classical polylogarithm Li_m(x) using series representation up to order n integer :: m integer :: m complex(kind=prec) :: x, res complex(kind=prec) :: x, res integer :: i integer(kind=ikin) :: i res=0._prec res=0._prec do i=1,PolylogInfinity do i=1,PolylogInfinity if(i**m.lt.0) return ! roll over if(i**m.lt.0) return ! roll over ... ...
 ... @@ -62,7 +62,7 @@ CONTAINS ... @@ -62,7 +62,7 @@ CONTAINS complex(kind=prec) :: x(:) complex(kind=prec) :: x(:) complex(kind=prec) :: res complex(kind=prec) :: res complex(kind=prec) :: t(size(x)) complex(kind=prec) :: t(size(x)) integer :: q, j, k integer(kind=ikin) :: q, j, k #ifdef MPL_CACHE #ifdef MPL_CACHE if (check_cache(m,x,res)) return if (check_cache(m,x,res)) return #endif #endif ... @@ -70,6 +70,9 @@ CONTAINS ... @@ -70,6 +70,9 @@ CONTAINS j = size(x) j = size(x) #ifdef DEBUG if(verb >= 70) print*, 'called MPL(',m,',',x,')' #endif if(size(m) /= size(x)) then if(size(m) /= size(x)) then print*, 'Error: m and x must have the same length' print*, 'Error: m and x must have the same length' end if end if ... @@ -79,13 +82,16 @@ CONTAINS ... @@ -79,13 +82,16 @@ CONTAINS do while(.true.) do while(.true.) res = t(1) res = t(1) q = q+1 q = q+1 if (q < 0) exit if(q**m(j) .eq. 0) exit t(j) = t(j) + x(j)**q/q**m(j) t(j) = t(j) + x(j)**q/q**m(j) do k=1,j-1 do k=1,j-1 if((k+q)**m(j-k) .eq. 0) exit t(j-k) = t(j-k) + t(j-k+1) * x(j-k)**(k+q)/(k+q)**m(j-k) t(j-k) = t(j-k) + t(j-k+1) * x(j-k)**(k+q)/(k+q)**m(j-k) enddo enddo if (mod(q,2) .eq. 1) then if (mod(q,2_ikin) .eq. 1) then if (abs(t(1)-res).lt.MPLdel) exit if (abs(t(1)-res).lt.MPLdel/100.) exit endif endif enddo enddo res = t(1) res = t(1) ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!