Bugfix: preserve sign of zero (#27)

---------

Co-authored-by: Alexander Voigt <alexander.voigt@physik.rwth-aachen.de>

Author: Expander

src/Li.jl

@@ -32,7 +32,7 @@ _reli(n::Integer, x::Float32) = oftype(x, _reli(n, Float64(x)))
 function _reli(n::Integer, x::Real)
     isnan(x) && return oftype(x, NaN)
     isinf(x) && return oftype(x, -Inf)
-    iszero(x) && return zero(x)
+    iszero(x) && return x
     x == one(x) && return zeta(n, typeof(x))
     x == -one(x) && return neg_eta(n, typeof(x))
 
@@ -120,7 +120,7 @@ function _li(n::Integer, z::Complex{T})::Complex{T} where T
 
     if iszero(imag(z))
         if real(z) <= one(typeof(real(z))) || n <= 0
-            return Complex(reli(n, real(z)))
+            return Complex(reli(n, real(z)), imag(z))
         else
             return Complex(reli(n, real(z)), -convert(T, pi)*inv_fac(n - 1, T)*log(real(z))^(n - 1))
         end

src/Li0.jl

@@ -18,4 +18,10 @@ julia> li0(BigFloat("0.25"))
 0.3333333333333333333333333333333333333333333333333333333333333333333333333333348
 ```
 """
-li0(z::Number) = z/(1 - z)
+function li0(z::Number)
+    if iszero(z)
+        z
+    else
+        z/(1 - z)
+    end
+end

src/Li1.jl

@@ -19,7 +19,9 @@ julia> reli1(BigFloat("0.5"))
 ```
 """
 function reli1(x::Real)
-    if x < one(x)
+    if x == zero(x)
+        x
+    elseif x < one(x)
         -log(one(x) - x)
     elseif x == one(x)
         oftype(x, Inf)
@@ -55,6 +57,12 @@ julia> li1(BigFloat("1.0") + 1im)
 -0.0 + 1.570796326794896619231321691639751442098584699687552910487472296153908203143099im
 ```
 """
-li1(z::Complex) = -clog1p(-z)
+function li1(z::Complex)
+    if iszero(z)
+        z
+    else
+        -clog1p(-z)
+    end
+end
 
 li1(z::Real) = li1(Complex(z))

src/Li2.jl

@@ -187,7 +187,7 @@ function _reli2(x::T)::T where T
     elseif x < zero(x)
         -reli2_approx(x/(x - one(x))) - one(x)/2*log1p(-x)^2
     elseif iszero(x)
-        zero(x)
+        x
     elseif x < one(x)/2
         reli2_approx(x)
     elseif x == one(x)/2
@@ -242,9 +242,9 @@ function _li2(z::Complex{T})::Complex{T} where T
 
     if iszero(iz)
         if rz <= one(T)
-            complex(reli2(rz))
+            Complex(reli2(rz), iz)
         else # Re(z) > 1
-            complex(reli2(rz), -convert(T, pi)*log(rz))
+            Complex(reli2(rz), -convert(T, pi)*log(rz))
         end
     else
         nz = abs2(z)
@@ -275,9 +275,9 @@ function _li2(z::Complex{BigFloat})::Complex{BigFloat}
 
     if iszero(iz)
         if rz <= one(BigFloat)
-            complex(reli2(rz))
+            Complex(reli2(rz), iz)
         else # Re(z) > 1
-            complex(reli2(rz), -convert(BigFloat, pi)*log(rz))
+            Complex(reli2(rz), -convert(BigFloat, pi)*log(rz))
         end
     else
         nz = abs2(z)

src/Li3.jl

@@ -84,7 +84,7 @@ function _reli3(x::Float64)::Float64
     elseif x < 0.0
         reli3_neg(x)
     elseif iszero(x)
-        0.0
+        x
     elseif x < 0.5
         reli3_pos(x)
     elseif x == 0.5
@@ -137,9 +137,9 @@ _li3(z::ComplexF32) = oftype(z, _li3(ComplexF64(z)))
 function _li3(z::ComplexF64)::ComplexF64
     if iszero(imag(z))
         if real(z) <= 1.0
-            reli3(real(z)) + 0.0im
+            Complex(reli3(real(z)), imag(z))
         else
-            reli3(real(z)) - pi*log(real(z))^2*0.5im
+            Complex(reli3(real(z)), -0.5*pi*log(real(z))^2)
         end
     else # Im(z) != 0
         nz::Float64  = abs(z)

src/Li4.jl

@@ -128,6 +128,8 @@ function _reli4(x::Float64)::Float64
 
     if x < 0.0
         rest + sgn*reli4_neg(x)
+    elseif iszero(x)
+        x
     elseif x < 0.5
         rest + sgn*reli4_half(x)
     elseif x < 0.8
@@ -172,9 +174,9 @@ _li4(z::ComplexF32) = oftype(z, _li4(ComplexF64(z)))
 function _li4(z::ComplexF64)::ComplexF64
     if iszero(imag(z))
         if real(z) <= 1.0
-            reli4(real(z)) + 0.0im
+            Complex(reli4(real(z)), imag(z))
         else
-            reli4(real(z)) - pi/6*log(real(z))^3*1.0im
+            Complex(reli4(real(z)), -pi/6*log(real(z))^3)
         end
     else
         nz::Float64  = abs(z)

src/Li5.jl

@@ -17,20 +17,22 @@ julia> li5(1.0 + 1.0im)
 """
 li5(z::Complex) = _li5(float(z))
 
+li5(z::Real) = li5(Complex(z))
+
 _li5(z::ComplexF16) = oftype(z, _li5(ComplexF32(z)))
 
 _li5(z::ComplexF32) = oftype(z, _li5(ComplexF64(z)))
 
 function _li5(z::ComplexF64)::ComplexF64
     if iszero(imag(z))
         if iszero(real(z))
-            return 0.0 + 0.0im
+            return z
         end
         if real(z) == 1.0
-            return ZETA5_F64 + 0.0im
+            return complex(ZETA5_F64)
         end
         if real(z) == -1.0
-            return -15/16*ZETA5_F64 + 0.0im
+            return complex(-15/16*ZETA5_F64)
         end
     end
 

src/Li6.jl

@@ -17,20 +17,22 @@ julia> li6(1.0 + 1.0im)
 """
 li6(z::Complex) = _li6(float(z))
 
+li6(z::Real) = li6(Complex(z))
+
 _li6(z::ComplexF16) = oftype(z, _li6(ComplexF32(z)))
 
 _li6(z::ComplexF32) = oftype(z, _li6(ComplexF64(z)))
 
 function _li6(z::ComplexF64)::ComplexF64
     if iszero(imag(z))
         if iszero(real(z))
-            return 0.0 + 0.0im
+            return z
         end
         if real(z) == 1.0
-            return ZETA6_F64 + 0.0im
+            return complex(ZETA6_F64)
         end
         if real(z) == -1.0
-            return -31/32*ZETA6_F64 + 0.0im
+            return complex(-31/32*ZETA6_F64)
         end
     end
 

test/Li.jl

@@ -55,6 +55,22 @@ end
             end
         end
 
+        # test signbit for 0.0 and -0.0 arguments
+        for T in (Float16, Float32, Float64)
+            @test  signbit(real(PolyLog.li(n, T(-0.0))))
+            @test !signbit(imag(PolyLog.li(n, T(-0.0))))
+            @test !signbit(real(PolyLog.li(n, T( 0.0))))
+            @test !signbit(imag(PolyLog.li(n, T( 0.0))))
+            @test !signbit(real(PolyLog.li(n, Complex{T}(0.0, 0.0))))
+            @test !signbit(imag(PolyLog.li(n, Complex{T}(0.0, 0.0))))
+            @test  signbit(real(PolyLog.li(n, Complex{T}(-0.0, 0.0))))
+            @test !signbit(imag(PolyLog.li(n, Complex{T}(-0.0, 0.0))))
+            @test !signbit(real(PolyLog.li(n, Complex{T}(0.0, -0.0))))
+            @test  signbit(imag(PolyLog.li(n, Complex{T}(0.0, -0.0))))
+            @test  signbit(real(PolyLog.li(n, Complex{T}(-0.0, -0.0))))
+            @test  signbit(imag(PolyLog.li(n, Complex{T}(-0.0, -0.0))))
+        end
+
         zeta = PolyLog.zeta(n, Float64)
 
         @test PolyLog.reli(n, 1.0) == zeta

test/Li0.jl

@@ -16,6 +16,22 @@
     test_function_on_data(PolyLog.li0, map(ComplexF16, cmpl_data), 1e-2, 1e-2)
     test_function_on_data(PolyLog.li0, map(Float16   , real_data), 1e-2, 1e-2)
 
+    # test signbit for 0.0 and -0.0 arguments
+    for T in (Float16, Float32, Float64, BigFloat)
+        @test  signbit(real(PolyLog.li0(T(-0.0))))
+        @test !signbit(imag(PolyLog.li0(T(-0.0))))
+        @test !signbit(real(PolyLog.li0(T( 0.0))))
+        @test !signbit(imag(PolyLog.li0(T( 0.0))))
+        @test !signbit(real(PolyLog.li0(Complex{T}(0.0, 0.0))))
+        @test !signbit(imag(PolyLog.li0(Complex{T}(0.0, 0.0))))
+        @test  signbit(real(PolyLog.li0(Complex{T}(-0.0, 0.0))))
+        @test !signbit(imag(PolyLog.li0(Complex{T}(-0.0, 0.0))))
+        @test !signbit(real(PolyLog.li0(Complex{T}(0.0, -0.0))))
+        @test  signbit(imag(PolyLog.li0(Complex{T}(0.0, -0.0))))
+        @test  signbit(real(PolyLog.li0(Complex{T}(-0.0, -0.0))))
+        @test  signbit(imag(PolyLog.li0(Complex{T}(-0.0, -0.0))))
+    end
+
     @test PolyLog.li0(2.0) ≈ -2.0
     @test PolyLog.li0(2.0f0) ≈ -2.0f0
     @test PolyLog.li0(Float16(2.0)) ≈ Float16(-2.0)