diff --git a/src/Constants.jl b/src/Constants.jl index 025c7a8..00e5711 100644 --- a/src/Constants.jl +++ b/src/Constants.jl @@ -1,4 +1,5 @@ -const zeta2F32 = 1.644934f0 # zeta(2) with Float32 precision +const zeta2F16 = Float16(1.64493) # zeta(2) with Float16 precision +const zeta2F32 = 1.644934f0 # zeta(2) with Float32 precision const zeta2 = 1.6449340668482264 const zeta3 = 1.2020569031595943 const zeta4 = 1.0823232337111382 diff --git a/src/Li2.jl b/src/Li2.jl index d8b8ae6..41c5e4a 100644 --- a/src/Li2.jl +++ b/src/Li2.jl @@ -235,70 +235,35 @@ li2(z::Real) = li2(Complex(z)) _li2(z::ComplexF16) = oftype(z, _li2(ComplexF32(z))) -function _li2(z::ComplexF32)::ComplexF32 - clog(z) = log(abs(z)) + angle(z)*1.0f0im +function _li2(z::Complex{T})::Complex{T} where T + clog(z) = T == Float64 ? complex(0.5*log(abs2(z)), angle(z)) : complex(log(abs(z)), angle(z)) rz, iz = reim(z) - if iz == 0.0f0 - if rz <= 1.0f0 - reli2(rz) + if iszero(iz) + if rz <= one(T) + complex(reli2(rz)) else # Re(z) > 1 - reli2(rz) - pi*clog(rz)*1.0f0im + complex(reli2(rz), -convert(T, pi)*log(rz)) end else nz = abs2(z) - if nz < eps(Float32) - z*(1.0f0 + 0.25f0*z) + if nz < eps(T) + z*(one(T) + one(T)/4*z) else - if rz <= 0.5f0 - if nz > 1.0f0 - -li2_approx(-clog(1.0f0 - inv(z))) - 0.5f0*clog(-z)^2 - zeta2F32 + if rz <= one(T)/2 + if nz > one(T) + -li2_approx(-clog(one(T) - inv(z))) - one(T)/2*clog(-z)^2 - zeta_2(T) else # |z|^2 <= 1 - li2_approx(-clog(1.0f0 - z)) + li2_approx(-clog(one(T) - z)) end else # Re(z) > 1/2 - if nz <= 2.0f0*rz + if nz <= 2*rz l = -clog(z) - -li2_approx(l) + l*clog(1.0f0 - z) + zeta2F32 + -li2_approx(l) + l*clog(one(T) - z) + zeta_2(T) else # |z|^2 > 2*Re(z) - -li2_approx(-clog(1.0f0 - inv(z))) - 0.5f0*clog(-z)^2 - zeta2F32 - end - end - end - end -end - -function _li2(z::ComplexF64)::ComplexF64 - clog(z) = 0.5*log(abs2(z)) + angle(z)*1.0im - - rz, iz = reim(z) - - if iz == 0.0 - if rz <= 1.0 - reli2(rz) - else # Re(z) > 1 - reli2(rz) - pi*log(rz)*1.0im - end - else - nz = abs2(z) - - if nz < eps(Float64) - z*(1.0 + 0.25*z) - else - if rz <= 0.5 - if nz > 1.0 - -li2_approx(-clog(1.0 - inv(z))) - 0.5*clog(-z)^2 - zeta2 - else # |z|^2 <= 1 - li2_approx(-clog(1.0 - z)) - end - else # Re(z) > 1/2 - if nz <= 2.0*rz - l = -clog(z) - -li2_approx(l) + l*clog(1.0 - z) + zeta2 - else # |z|^2 > 2*Re(z) - -li2_approx(-clog(1.0 - inv(z))) - 0.5*clog(-z)^2 - zeta2 + -li2_approx(-clog(one(T) - inv(z))) - one(T)/2*clog(-z)^2 - zeta_2(T) end end end diff --git a/src/Zeta.jl b/src/Zeta.jl index adccae4..d1a70ed 100644 --- a/src/Zeta.jl +++ b/src/Zeta.jl @@ -99,3 +99,16 @@ function zeta_big(n::Integer)::BigFloat ccall((:mpfr_zeta, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[]) return z end + +# returns pre-computed value zeta(2) for given type T +function zeta_2(::Type{T})::T where T + if T == Float16 + zeta2F16 + elseif T == Float32 + zeta2F32 + elseif T == Float64 + zeta2 + else + zeta(2, T) + end +end diff --git a/test/Zeta.jl b/test/Zeta.jl index 52fe4d6..d6cf37c 100644 --- a/test/Zeta.jl +++ b/test/Zeta.jl @@ -39,3 +39,12 @@ @test PolyLog.zeta( 8, BigFloat) ≈ big(pi)^8/9450 rtol=10*eps(BigFloat) @test PolyLog.zeta( 10, BigFloat) ≈ big(pi)^10/93555 rtol=10*eps(BigFloat) end + +@testset "zeta_2" begin + zeta2F64 = 1.6449340668482264 + + @test PolyLog.zeta_2(Float16) == Float16(zeta2F64) + @test PolyLog.zeta_2(Float32) == Float32(zeta2F64) + @test PolyLog.zeta_2(Float64) == zeta2F64 + @test PolyLog.zeta_2(BigFloat) == PolyLog.zeta(2, BigFloat) +end