...
 
Commits (2)
......@@ -400,6 +400,19 @@ CONTAINS
END FUNCTION
FUNCTION MYLOG(x)
complex(kind=prec) :: x, mylog
if (abs(aimag(x)) < zero) then
if (real(x) > 0) then
mylog = cmplx(log(real(+x,kind=prec)), kind=prec)
else
mylog = cmplx(log(real(-x,kind=prec)), pi, kind=prec)
endif
else
mylog = log(x)
endif
END FUNCTION
RECURSIVE FUNCTION polylog1(m,x) result(res)
! computes the polylog
......@@ -433,7 +446,7 @@ CONTAINS
else if (abs(x) .gt. 1) then
inv = 1._prec/x
res = (-1)**(m-1)*polylog(m,inv) &
- (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*log(-x)/2._prec/pi) / factorial(m)
- (2._prec*pi*i_)**m * bernoulli_polynomial(m, 0.5_prec-i_*mylog(-x)/2._prec/pi) / factorial(m)
else if(m == 2) then
res = dilog(x)
else if(m == 3) then
......@@ -489,7 +502,7 @@ CONTAINS
plog1 = plog1 - b%i0*i_*pi
endif
else
plog1 = log(1.-a%c/b%c)
plog1 = mylog(1.-a%c/b%c)
endif
END FUNCTION
......