Commit 3215f5fb authored by ulrich_y's avatar ulrich_y

Removed old ieps framework

parent 61308a0a
......@@ -4,7 +4,6 @@ MODULE globals
integer, parameter :: prec = selected_real_kind(15,32)
integer, parameter :: GPLInfinity = 30 ! the default outermost expansion order for MPLs
real, parameter :: epsilon = 1e-20 ! used for the small imaginary part
real, parameter :: zero = 1e-15 ! values smaller than this count as zero
real, parameter :: pi = 3.14159265358979323846
......
......@@ -124,19 +124,15 @@ CONTAINS
END FUNCTION Li2
RECURSIVE FUNCTION dilog(x) result(res)
! evaluates dilog for any argument
! evaluates dilog for any argument |x|<1
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(x) <= 1.0) then
if(abs(aimag(x)) < zero ) then
res = Li2(real(x))
else
res = naive_polylog(2,x)
endif
if(abs(aimag(x)) < zero ) then
res = Li2(real(x))
else
res = -dilog(1/x) - (pi**2) /6 - log(add_ieps(-x))**2 / 2
end if
res = naive_polylog(2,x)
endif
END FUNCTION dilog
FUNCTION Li3(x)
......@@ -244,18 +240,14 @@ CONTAINS
END FUNCTION Li3
FUNCTION trilog(x) result(res)
! evaluates trilog for any argument
! evaluates trilog for any argument |x|<1
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(x) <= 1.0) then
if(abs(aimag(x)) < zero ) then
res = Li3(real(x))
else
res = naive_polylog(3,x)
endif
if(abs(aimag(x)) < zero ) then
res = Li3(real(x))
else
res = naive_polylog(3,sub_ieps(x)**(-1)) - (log(-sub_ieps(x)))**3/6 - pi**2/6 * log(-sub_ieps(x))
end if
res = naive_polylog(3,x)
endif
END FUNCTION trilog
FUNCTION BERNOULLI_POLYNOMIAL(n, x) result(res)
......
......@@ -126,18 +126,6 @@ CONTAINS
res = merge(1,n*factorial(n-1),n==0)
END FUNCTION factorial
FUNCTION add_ieps(x) result(res)
! adds small imaginary part to x
complex(kind=prec) :: x, res
res = x + (0.0,epsilon)
END FUNCTION add_ieps
FUNCTION sub_ieps(x) result(res)
! subtracts small imaginary part to x
complex(kind=prec) :: x, res
res = x - (0.0,epsilon)
END FUNCTION sub_ieps
FUNCTION shuffle_with_zero(a) result(res)
! rows of result are shuffles of a with 0
type(inum) :: a(:)
......
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