preserve sign of 0 for real-valued arguments

Expander · Feb 21, 2024 · cfea54c · cfea54c
1 parent fad19d6
commit cfea54c
Show file tree

Hide file tree

Showing 15 changed files with 101 additions and 15 deletions.
diff --git a/src/Li.jl b/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))
 

diff --git a/src/Li1.jl b/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))
diff --git a/src/Li2.jl b/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

diff --git a/src/Li3.jl b/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)))
         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)

diff --git a/src/Li4.jl b/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)))
         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)

diff --git a/src/Li5.jl b/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
 

diff --git a/src/Li6.jl b/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
 

diff --git a/test/Li.jl b/test/Li.jl
@@ -55,6 +55,14 @@ 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))))
+        end
+
         zeta = PolyLog.zeta(n, Float64)
 
         @test PolyLog.reli(n, 1.0) == zeta

diff --git a/test/Li0.jl b/test/Li0.jl
@@ -16,6 +16,14 @@
     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))))
+    end
+
     @test PolyLog.li0(2.0) ≈ -2.0
     @test PolyLog.li0(2.0f0) ≈ -2.0f0
     @test PolyLog.li0(Float16(2.0)) ≈ Float16(-2.0)

diff --git a/test/Li1.jl b/test/Li1.jl
@@ -17,6 +17,16 @@
     test_function_on_data(PolyLog.li1  , map(ComplexF16, cmpl_data), 1e-2, 1e-2)
     test_function_on_data(PolyLog.reli1, 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(PolyLog.reli1(T(-0.0)))
+        @test !signbit(PolyLog.reli1(T( 0.0)))
+        @test  signbit(real(PolyLog.li1(T(-0.0))))
+        @test !signbit(imag(PolyLog.li1(T(-0.0))))
+        @test !signbit(real(PolyLog.li1(T( 0.0))))
+        @test !signbit(imag(PolyLog.li1(T( 0.0))))
+    end
+
     @test PolyLog.reli1(-1.0) ≈ -log(2.0)
     @test PolyLog.reli1(-1.0f0) ≈ -log(2.0f0)
     @test PolyLog.reli1(Float16(-1.0)) ≈ Float16(-log(2.0))

diff --git a/test/Li2.jl b/test/Li2.jl
@@ -18,6 +18,16 @@
     test_function_on_data(PolyLog.li2  , filter_ComplexF16(map(ComplexF16, cmpl_data)), 1e-2, 1e-2)
     test_function_on_data(PolyLog.reli2, 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(PolyLog.reli2(T(-0.0)))
+        @test !signbit(PolyLog.reli2(T( 0.0)))
+        @test  signbit(real(PolyLog.li2(T(-0.0))))
+        @test !signbit(imag(PolyLog.li2(T(-0.0))))
+        @test !signbit(real(PolyLog.li2(T( 0.0))))
+        @test !signbit(imag(PolyLog.li2(T( 0.0))))
+    end
+
     zeta2 = 1.6449340668482264
 
     @test PolyLog.reli2(1.0) == zeta2

diff --git a/test/Li3.jl b/test/Li3.jl
@@ -11,6 +11,16 @@
     test_function_on_data(PolyLog.li3  , filter_ComplexF16(map(ComplexF16, cmpl_data)), 1e-2, 1e-2)
     test_function_on_data(PolyLog.reli3, map(Float16, real_data), 1e-2, 1e-2)
 
+    # test signbit for 0.0 and -0.0 arguments
+    for T in (Float16, Float32, Float64)
+        @test  signbit(PolyLog.reli3(T(-0.0)))
+        @test !signbit(PolyLog.reli3(T( 0.0)))
+        @test  signbit(real(PolyLog.li3(T(-0.0))))
+        @test !signbit(imag(PolyLog.li3(T(-0.0))))
+        @test !signbit(real(PolyLog.li3(T( 0.0))))
+        @test !signbit(imag(PolyLog.li3(T( 0.0))))
+    end
+
     zeta3 = 1.2020569031595943
 
     @test PolyLog.reli3(1.0) == zeta3

diff --git a/test/Li4.jl b/test/Li4.jl
@@ -11,6 +11,16 @@
     test_function_on_data(PolyLog.li4  , map(ComplexF16, cmpl_data), 1e-2, 1e-2)
     test_function_on_data(PolyLog.reli4, map(Float16   , real_data), 1e-2, 1e-2)
 
+    # test signbit for 0.0 and -0.0 arguments
+    for T in (Float16, Float32, Float64)
+        @test  signbit(PolyLog.reli4(T(-0.0)))
+        @test !signbit(PolyLog.reli4(T( 0.0)))
+        @test  signbit(real(PolyLog.li4(T(-0.0))))
+        @test !signbit(imag(PolyLog.li4(T(-0.0))))
+        @test !signbit(real(PolyLog.li4(T( 0.0))))
+        @test !signbit(imag(PolyLog.li4(T( 0.0))))
+    end
+
     zeta4 = 1.0823232337111382
 
     @test PolyLog.reli4(1.0) == zeta4

diff --git a/test/Li5.jl b/test/Li5.jl
@@ -5,6 +5,14 @@
     test_function_on_data(PolyLog.li5, map(ComplexF32, cmpl_data), 1e-6, 1e-6)
     test_function_on_data(PolyLog.li5, map(ComplexF16, cmpl_data), 1e-2, 1e-2)
 
+    # test signbit for 0.0 and -0.0 arguments
+    for T in (Float16, Float32, Float64)
+        @test  signbit(real(PolyLog.li5(T(-0.0))))
+        @test !signbit(imag(PolyLog.li5(T(-0.0))))
+        @test !signbit(real(PolyLog.li5(T( 0.0))))
+        @test !signbit(imag(PolyLog.li5(T( 0.0))))
+    end
+
     zeta5 = 1.0369277551433699
 
     @test PolyLog.li5(1.0 + 0.0im) == zeta5

diff --git a/test/Li6.jl b/test/Li6.jl
@@ -5,6 +5,14 @@
     test_function_on_data(PolyLog.li6, map(ComplexF32, cmpl_data), 1e-6, 1e-6)
     test_function_on_data(PolyLog.li6, map(ComplexF16, cmpl_data), 1e-2, 1e-2)
 
+    # test signbit for 0.0 and -0.0 arguments
+    for T in (Float16, Float32, Float64)
+        @test  signbit(real(PolyLog.li6(T(-0.0))))
+        @test !signbit(imag(PolyLog.li6(T(-0.0))))
+        @test !signbit(real(PolyLog.li6(T( 0.0))))
+        @test !signbit(imag(PolyLog.li6(T( 0.0))))
+    end
+
     zeta6 = 1.0173430619844491
 
     @test PolyLog.li6(1.0 + 0.0im) == zeta6