add complex log1p (fix JuliaLang#3141)

Sacha0 · May 21, 2014 · 2366342 · 2366342
1 parent 2bc0003
commit 2366342
Show file tree

Hide file tree

Showing 4 changed files with 45 additions and 4 deletions.
diff --git a/NEWS.md b/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]).

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

diff --git a/base/math.jl b/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
 

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