Merge branch 'main' into bug-imprecise-complex-li2-for-small-z

Expander · Feb 20, 2024 · cb8de2c · cb8de2c
2 parents 24bc618 + 7b03b0b
commit cb8de2c
Show file tree

Hide file tree

Showing 7 changed files with 21 additions and 15 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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)