From 09592838d02d3fe6257aaba011ef9900130f02aa Mon Sep 17 00:00:00 2001 From: Yannick Ulrich Date: Thu, 3 Oct 2019 21:08:18 +0200 Subject: [PATCH] Fixed branch cut of polylog again This fixes the problem of 3b2d170d and bdc545ff for good. In the previous version I had a conjg in there, but that turned out to be nonsense. To see why, see the comment in 3b2d170d (spoiler alert: it has to do with a peculiarity of IEEE 754) --- src/maths_functions.f90 | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/maths_functions.f90 b/src/maths_functions.f90 index 0bbb1ec..3ff3246 100644 --- a/src/maths_functions.f90 +++ b/src/maths_functions.f90 @@ -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 -- GitLab