Commit f652eea6 authored by ulrich_y's avatar ulrich_y

polylog2 is now careful with the branch cut

parent 671a9a2c
......@@ -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
......
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