add complex log1p (fix #3141)

Author: stevengj

NEWS.md

@@ -144,6 +144,8 @@ Library improvements
 
     * `rand` now supports arbitrary `Ranges` arguments ([#5059]).
 
+    * `expm1` and `log1p` now support complex arguments ([#3141]).
+
   * `String` improvements
 
     * Triple-quoted regex strings, `r"""..."""` ([#4934]).

base/complex.jl

@@ -329,7 +329,7 @@ log2(z::Complex) = log(z)/oftype(real(z),0.6931471805599453)
 function exp(z::Complex)
     zr, zi = reim(z)
     if isnan(zr)
-        Complex(zr, zi==zero(zi) ? zi : zr)
+        Complex(zr, zi==0 ? zi : zr)
     elseif !isfinite(zi)
         if zr == inf(zr)
             Complex(-zr, nan(zr))
@@ -351,7 +351,7 @@ end
 function expm1(z::Complex)
     zr,zi = reim(z)
     if isnan(zr)
-        Complex(zr, zi==zero(zi) ? zi : zr)
+        Complex(zr, zi==0 ? zi : zr)
     elseif !isfinite(zi)
         if zr == inf(zr)
             Complex(-zr, nan(zr))
@@ -362,7 +362,7 @@ function expm1(z::Complex)
         end
     else
         erm1 = expm1(zr)        
-        if zi == zero(zi)
+        if zi == 0
             Complex(erm1, zi)
         else            
             er = erm1+one(erm1)
@@ -372,6 +372,24 @@ function expm1(z::Complex)
     end
 end
 
+function log1p{T}(z::Complex{T})
+    zr,zi = reim(z)
+    if isfinite(zr)
+        isinf(zi) && return log(z)
+        # This is based on a well-known trick for log1p of real z,
+        # allegedly due to Kahan, only modified to handle real(u) <= 0
+        # differently to avoid inaccuracy near z==-2 and for correct branch cut
+        u = float(one(T)) + z
+        u == 1 ? convert(typeof(u), z) : real(u) <= 0 ? log(u) : log(u)*z/(u-1)
+    elseif isnan(zr)
+        Complex(zr, zr)
+    elseif isfinite(zi)
+        Complex(inf(T), copysign(zr > 0 ? zero(T) : convert(T, pi), zi))
+    else
+        Complex(inf(T), nan(T))
+    end
+end
+
 function ^{T<:FloatingPoint}(z::Complex{T}, p::Complex{T})
     if p==2 #square
         zr, zi = reim(z)

base/math.jl

@@ -21,7 +21,7 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
 import Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
              acos, atan, asinh, acosh, atanh, sqrt, log2, log10,
              max, min, minmax, ceil, floor, trunc, round, ^, exp2,
-             exp10, expm1
+             exp10, expm1, log1p
 
 import Core.Intrinsics: nan_dom_err, sqrt_llvm, box, unbox, powi_llvm
 

test/complex.jl

@@ -151,6 +151,27 @@
 @test_approx_eq expm1(complex(-10.0, 10.0)) exp(complex(-10.0, 10.0))-1
 @test_approx_eq expm1(complex(-10.0,-10.0)) exp(complex(-10.0,-10.0))-1
 
+# log1p:
+@test isequal(log1p(complex(Inf, 3)), complex(Inf, +0.))
+@test isequal(log1p(complex(Inf, -3)), complex(Inf, -0.))
+@test isequal(log1p(complex(-Inf, 3)), complex(Inf, +pi))
+@test isequal(log1p(complex(-Inf, -3)), complex(Inf, -pi))
+@test isequal(log1p(complex(Inf, NaN)), complex(Inf, NaN))
+@test isequal(log1p(complex(NaN, 0)), complex(NaN, NaN))
+@test isequal(log1p(complex(0, NaN)), complex(NaN, NaN))
+@test isequal(log1p(complex(-1, +0.)), complex(-Inf, +0.))
+@test isequal(log1p(complex(-1, -0.)), complex(-Inf, -0.))
+@test isequal(log1p(complex(-2, 1e-10)), log(1 + complex(-2, 1e-10)))
+@test isequal(log1p(complex(1, Inf)), complex(Inf, pi/2))
+@test isequal(log1p(complex(1, -Inf)), complex(Inf, -pi/2))
+import Base.Math.@horner
+for z in (1e-10+1e-9im, 1e-10-1e-9im, -1e-10+1e-9im, -1e-10-1e-9im)
+    @test_approx_eq log1p(z) @horner(z, 0, 1, -0.5, 1/3, -0.25, 0.2)
+end
+for z in (15+4im, 0.2+3im, 0.08-0.9im)
+    @test_approx_eq log1p(z) log(1+z)
+end
+
 
 # ^ (cpow)
 #  equivalent to exp(y*log(x))