diff --git a/src/maths_functions.f90 b/src/maths_functions.f90 index 04cae629f56b187891093db13b203cd044a8010e..cf7180b7f7df4bd155143549c13057f785b8ef34 100644 --- a/src/maths_functions.f90 +++ b/src/maths_functions.f90 @@ -14,7 +14,7 @@ MODULE maths_functions 1.0083492773819228268397975498497967595998635605652_prec, 1.0040773561979443393786852385086524652589607906499_prec, & 1.0020083928260822144178527692324120604856058513949_prec, 1.0009945751278180853371459589003190170060195315645_prec /) type el - type(inum) :: c + complex(kind=prec) :: c complex(kind=prec) ans end type el @@ -398,38 +398,38 @@ CONTAINS ! computes the polylog integer :: m - type(inum) :: x, inv + complex(kind=prec) :: x, inv complex(kind=prec) :: res integer i #ifdef DEBUG - if(verb >= 70) print*, 'called polylog(',m,',',x%c,x%i0,')' + if(verb >= 70) print*, 'called polylog(',m,',',x,')' #endif #ifndef NOCACHE if (m.le.5) then do i=1,plcachesize(m) - if( abs(cache(m,i)%c%c-x%c).lt.zero ) then + if( abs(cache(m,i)%c-x).lt.zero ) then res = cache(m,i)%ans return endif enddo endif #endif - if ((m.le.9).and.(abs(x%c-1.).lt.zero)) then + if ((m.le.9).and.(abs(x-1.).lt.zero)) then res = zeta(m) - else if ((m.le.9).and.(abs(x%c+1._prec).lt.zero)) then + else if ((m.le.9).and.(abs(x+1._prec).lt.zero)) then res = -(1._prec - 2._prec**(1-m))*zeta(m) else if (abs(x) .gt. 1) then - inv = inum(1._prec/x%c, x%i0) + inv = 1._prec/x res = (-1)**(m-1)*polylog(m,inv) & - - (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*conjg(log(-x%c))/2._prec/pi) / factorial(m) + - (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*conjg(log(-x))/2._prec/pi) / factorial(m) else if(m == 2) then - res = dilog(x%c) + res = dilog(x) else if(m == 3) then - res = trilog(x%c) + res = trilog(x) else - res = naive_polylog(m,x%c) + res = naive_polylog(m,x) end if #ifndef NOCACHE @@ -449,8 +449,17 @@ CONTAINS type(inum) :: x, y integer m complex(kind=prec) :: res - !TODO!! - res=polylog1(m,inum(x%c/y%c,di0)) + res=polylog1(m,x%c/y%c) + if ( (abs(aimag(x)).lt.zero).and.(abs(aimag(y)).lt.zero) ) then + ! Both arguments are real, only here does the ieps matter + ! FIXME this is rather ugly.. + if (real(x).gt.real(y) .and. real(y).gt. 0) then + ! Force the sign to be -b%i0 + res = cmplx(real(res), sign(aimag(res), -y%i0*1._prec), kind=prec) + elseif(real(x).lt.real(y) .and. real(y).lt. 0) then + res = cmplx(real(res), sign(aimag(res), +y%i0*1._prec), kind=prec) + endif + endif END FUNCTION POLYLOG2