Commit ef6e73d9 authored by ulrich_y's avatar ulrich_y

Prevent over- and under-flow in MPL

parent 41356ea9
......@@ -4,6 +4,9 @@ MODULE mpl_module
use ieps
implicit none
real(kind=prec), parameter :: underflowalert = 1.e-250_prec
real(kind=prec), parameter :: overflowalert = 1.e+250_prec
#ifdef MPL_CACHE
! Max weight, no. cache
integer,parameter :: cachesize(2) = (/ 4, 2500 /)
......@@ -61,8 +64,8 @@ CONTAINS
integer :: m(:)
complex(kind=prec) :: x(:)
complex(kind=prec) :: res
complex(kind=prec) :: t(size(x))
integer(kind=ikin) :: q, j, k
complex(kind=prec) :: t(size(x)), cpow
integer(kind=ikin) :: q, j, k, ipow
#ifdef MPL_CACHE
if (check_cache(m,x,res)) return
#endif
......@@ -83,11 +86,22 @@ CONTAINS
res = t(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)
cpow = x(j)**q
ipow = q**m(j)
if(ipow .eq. 0) exit
if (abs(cpow).lt.underflowalert) exit
if (abs(cpow).gt. overflowalert) exit
t(j) = t(j) + cpow/ipow
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)
ipow = (k+q)**m(j-k)
cpow = x(j-k)**(k+q)
if(ipow .eq. 0) exit
if (abs(cpow).lt.underflowalert) exit
if (abs(cpow).gt. overflowalert) exit
t(j-k) = t(j-k) + t(j-k+1) * cpow/ipow
enddo
if (mod(q,2_ikin) .eq. 1) then
......@@ -104,8 +118,8 @@ CONTAINS
return
endif
cacheind(j) = cacheind(j) + 1
cache(j,cacheind(j)) = el( reshape(m, (/cachesize(1)/), pad=[ 0 ]), &
reshape(x, (/cachesize(1)/), pad=[cmplx(0.)]), &
cache(j,cacheind(j)) = el( reshape(m, (/cachesize(1)/), pad=[ 0 ]), &
reshape(x, (/cachesize(1)/), pad=[cmplx(0.,kind=prec)]), &
res )
endif
#endif
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment