diff --git a/src/Li1.jl b/src/Li1.jl index 752a93e..1b6b697 100644 --- a/src/Li1.jl +++ b/src/Li1.jl @@ -55,6 +55,6 @@ julia> li1(BigFloat("1.0") + 1im) -0.0 + 1.570796326794896619231321691639751442098584699687552910487472296153908203143099im ``` """ -li1(z::Complex) = -clog(one(z) - z) +li1(z::Complex) = -clog1p(-z) li1(z::Real) = li1(Complex(z)) diff --git a/src/Li2.jl b/src/Li2.jl index ba5f2a0..16142b6 100644 --- a/src/Li2.jl +++ b/src/Li2.jl @@ -107,7 +107,7 @@ function li2_approx(z::Complex{T})::Complex{T} where T if abs2(z) < (99/100)^2 li2_approx_taylor(z) else - li2_approx_unity(-clog(one(T) - z)) + li2_approx_unity(-clog1p(-z)) end end @@ -254,16 +254,16 @@ function _li2(z::Complex{T})::Complex{T} where T else if rz <= one(T)/2 if nz > one(T) - -li2_approx(-clog(one(T) - inv(z))) - one(T)/2*clog(-z)^2 - zeta2(T) + -li2_approx(-clog1p(-inv(z))) - one(T)/2*clog(-z)^2 - zeta2(T) else # |z|^2 <= 1 - li2_approx(-clog(one(T) - z)) + li2_approx(-clog1p(-z)) end else # Re(z) > 1/2 if nz <= 2*rz l = -clog(z) - -li2_approx(l) + l*clog(one(T) - z) + zeta2(T) + -li2_approx(l) + l*clog1p(-z) + zeta2(T) else # |z|^2 > 2*Re(z) - -li2_approx(-clog(one(T) - inv(z))) - one(T)/2*clog(-z)^2 - zeta2(T) + -li2_approx(-clog1p(-inv(z))) - one(T)/2*clog(-z)^2 - zeta2(T) end end end @@ -293,7 +293,7 @@ function _li2(z::Complex{BigFloat})::Complex{BigFloat} end else # Re(z) > 1/2 if nz <= 2*rz - -li2_approx(one(BigFloat) - z) - clog(z)*clog(one(BigFloat) - z) + zeta2(BigFloat) + -li2_approx(one(BigFloat) - z) - clog(z)*clog1p(-z) + zeta2(BigFloat) else # |z|^2 > 2*Re(z) -li2_approx(inv(z)) - one(BigFloat)/2*clog(-z)^2 - zeta2(BigFloat) end diff --git a/src/Li3.jl b/src/Li3.jl index bc981d9..045666a 100644 --- a/src/Li3.jl +++ b/src/Li3.jl @@ -179,11 +179,11 @@ function _li3(z::ComplexF64)::ComplexF64 ) (u, rest) = if nz <= 1.0 - (-clog(1.0 - z), 0.0 + 0.0im) + (-clog1p(-z), 0.0 + 0.0im) else # |z|^2 > 1 arg = pz > 0.0 ? pz - pi : pz + pi lmz = lnz + arg*im # clog(z) - (-clog(1.0 - inv(z)), -lmz*(1/6*lmz*lmz + ZETA2_F64)) + (-clog1p(-inv(z)), -lmz*(1/6*lmz*lmz + ZETA2_F64)) end u2 = u*u diff --git a/src/Li4.jl b/src/Li4.jl index 7682d20..a8995a2 100644 --- a/src/Li4.jl +++ b/src/Li4.jl @@ -219,12 +219,12 @@ function _li4(z::ComplexF64)::ComplexF64 ) (u, rest, sgn) = if nz <= 1.0 - (-clog(1.0 - z), 0.0 + 0.0im, 1.0) + (-clog1p(-z), 0.0 + 0.0im, 1.0) else # |z|^2 > 1 arg = pz > 0.0 ? pz - pi : pz + pi lmz = lnz + arg*im # clog(z) lmz2 = lmz*lmz - (-clog(1.0 - inv(z)), -7/4*ZETA4_F64 + lmz2*(-0.5*ZETA2_F64 - 1/24*lmz2), -1.0) + (-clog1p(-inv(z)), -7/4*ZETA4_F64 + lmz2*(-0.5*ZETA2_F64 - 1/24*lmz2), -1.0) end u2 = u*u diff --git a/src/Li5.jl b/src/Li5.jl index a3de396..cd1c4bd 100644 --- a/src/Li5.jl +++ b/src/Li5.jl @@ -74,12 +74,12 @@ function _li5(z::ComplexF64)::ComplexF64 ) (u, rest) = if nz <= 1.0 - (-clog(1.0 - z), 0.0 + 0.0im) + (-clog1p(-z), 0.0 + 0.0im) else # nz > 1.0 arg = pz > 0.0 ? pz - pi : pz + pi lmz = lnz + arg*im # clog(z) lmz2 = lmz*lmz - (-clog(1.0 - inv(z)), lmz*(-7/4*ZETA4_F64 - lmz2*(1/6*ZETA2_F64 + 1/120*lmz2))) + (-clog1p(-inv(z)), lmz*(-7/4*ZETA4_F64 - lmz2*(1/6*ZETA2_F64 + 1/120*lmz2))) end u2 = u*u diff --git a/src/Li6.jl b/src/Li6.jl index 34a8d48..2fcc5ae 100644 --- a/src/Li6.jl +++ b/src/Li6.jl @@ -72,12 +72,12 @@ function _li6(z::ComplexF64)::ComplexF64 ) (u, rest, sgn) = if nz <= 1.0 - (-clog(1.0 - z), 0.0 + 0.0im, 1.0) + (-clog1p(-z), 0.0 + 0.0im, 1.0) else # nz > 1.0 arg = pz > 0.0 ? pz - pi : pz + pi lmz = lnz + arg*im # clog(z) lmz2 = lmz*lmz - (-clog(1.0 - inv(z)), -31/16*ZETA6_F64 + lmz2*(-7/8*ZETA4_F64 + lmz2*(-1/24*ZETA2_F64 - 1/720*lmz2)), -1.0) + (-clog1p(-inv(z)), -31/16*ZETA6_F64 + lmz2*(-7/8*ZETA4_F64 + lmz2*(-1/24*ZETA2_F64 - 1/720*lmz2)), -1.0) end u2 = u*u diff --git a/src/Log.jl b/src/Log.jl index adbb63a..2eeca1c 100644 --- a/src/Log.jl +++ b/src/Log.jl @@ -10,6 +10,12 @@ function clog(x::Real) log(x) end +# returns complex logarithm of 1+z; +# handles case when imag(z) == -0.0 +function clog1p(z::Complex) + clog(one(z) + z) +end + # returns |ln(x)|^2 for all x function ln_sqr(x::Real)::Real if x < zero(x)