Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adds an expm1 function for complex values. #6539

Merged
merged 2 commits into from
May 20, 2014
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 42 additions & 10 deletions base/complex.jl
Original file line number Diff line number Diff line change
Expand Up @@ -328,16 +328,48 @@ log2(z::Complex) = log(z)/oftype(real(z),0.6931471805599453)

function exp(z::Complex)
zr, zi = reim(z)
if isfinite(zr) && !isfinite(zi) return Complex(oftype(zr, NaN), oftype(zi, NaN)) end
if zr==Inf && zi==0 return Complex(zr, zi) end
if zr==-Inf && !isfinite(zi) return Complex(-zero(zr), copysign(zero(zi), zi)) end
if zr==Inf && !isfinite(zi) return Complex(-zr, oftype(zr, NaN)) end
if isnan(zr) return Complex(zr, zi==0 ? zi : zr) end
er = exp(zr)
zi==0 && return Complex(er, zi)
wr = er*(isfinite(zi) ? cos(zi) : zi)
wi = er*(isfinite(zi) ? sin(zi) : zi)
Complex(wr, wi)
if isnan(zr)
Complex(zr, zi==zero(zi) ? zi : zr)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not important but just FYI, zero() is not necessary in this context (since it's an argument and not a return value).

elseif !isfinite(zi)
if zr == inf(zr)
Complex(-zr, nan(zr))
elseif zr == -inf(zr)
Complex(-zero(zr), copysign(zero(zi), zi))
else
Complex(nan(zr), nan(zi))
end
else
er = exp(zr)
if zi == zero(zi)
Complex(er, zi)
else
Complex(er*cos(zi), er*sin(zi))
end
end
end

function expm1(z::Complex)
zr,zi = reim(z)
if isnan(zr)
Complex(zr, zi==zero(zi) ? zi : zr)
elseif !isfinite(zi)
if zr == inf(zr)
Complex(-zr, nan(zr))
elseif zr == -inf(zr)
Complex(-one(zr), copysign(zero(zi), zi))
else
Complex(nan(zr), nan(zi))
end
else
erm1 = expm1(zr)
if zi == zero(zi)
Complex(erm1, zi)
else
er = erm1+one(erm1)
wr = isfinite(er) ? erm1 - 2.0*er*(sin(0.5*zi))^2 : er*cos(zi)
Complex(wr, er*sin(zi))
end
end
end

function ^{T<:FloatingPoint}(z::Complex{T}, p::Complex{T})
Expand Down
3 changes: 2 additions & 1 deletion base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,8 @@ 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
max, min, minmax, ceil, floor, trunc, round, ^, exp2,
exp10, expm1

import Core.Intrinsics: nan_dom_err, sqrt_llvm, box, unbox, powi_llvm

Expand Down
49 changes: 49 additions & 0 deletions test/complex.jl
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,55 @@
@test isequal(exp(complex( NaN, 5.0)), complex( NaN, NaN))
@test isequal(exp(complex( NaN, NaN)), complex( NaN, NaN))

# expm1:
# expm1(conj(z)) = conj(expm1(z))

@test isequal(expm1(complex( 0.0, 0.0)), complex(0.0, 0.0))
@test isequal(expm1(complex( 0.0,-0.0)), complex(0.0,-0.0))
@test isequal(expm1(complex( 0.0, Inf)), complex(NaN, NaN))
@test isequal(expm1(complex( 0.0,-Inf)), complex(NaN, NaN))
@test isequal(expm1(complex( 0.0, NaN)), complex(NaN, NaN))

@test isequal(expm1(complex(-0.0, 0.0)), complex(-0.0, 0.0))
@test isequal(expm1(complex(-0.0,-0.0)), complex(-0.0,-0.0))

@test isequal(expm1(complex( 5.0, Inf)), complex(NaN, NaN))

@test isequal(expm1(complex( Inf, 0.0)), complex(Inf, 0.0))
@test isequal(expm1(complex( Inf,-0.0)), complex(Inf,-0.0))
@test isequal(expm1(complex( Inf, 5.0)), complex(cos(5.0)*Inf,sin(5.0)* Inf))
@test isequal(expm1(complex( Inf,-5.0)), complex(cos(5.0)*Inf,sin(5.0)*-Inf))
@test isequal(expm1(complex( Inf, NaN)), complex(-Inf, NaN))
@test isequal(expm1(complex( Inf, Inf)), complex(-Inf, NaN))
@test isequal(expm1(complex( Inf,-Inf)), complex(-Inf, NaN))

@test isequal(expm1(complex(-Inf, 0.0)), complex(-1.0, 0.0))
@test isequal(expm1(complex(-Inf,-0.0)), complex(-1.0,-0.0))
@test isequal(expm1(complex(-Inf, 5.0)), complex(-1.0,sin(5.0)* 0.0))
@test isequal(expm1(complex(-Inf,-5.0)), complex(-1.0,sin(5.0)*-0.0))
@test isequal(expm1(complex(-Inf, Inf)), complex(-1.0, 0.0))
@test isequal(expm1(complex(-Inf,-Inf)), complex(-1.0,-0.0))
@test isequal(expm1(complex(-Inf, NaN)), complex(-1.0, 0.0))

@test isequal(expm1(complex( NaN, 0.0)), complex( NaN, 0.0))
@test isequal(expm1(complex( NaN,-0.0)), complex( NaN,-0.0))
@test isequal(expm1(complex( NaN, 5.0)), complex( NaN, NaN))
@test isequal(expm1(complex( NaN, NaN)), complex( NaN, NaN))

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These tests are great but they don't seem to test any "normal" argument values.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, I've added some more tests.

@test isequal(expm1(complex( 1e-20, 0.0)), complex(expm1( 1e-20), 0.0))
@test isequal(expm1(complex(-1e-20, 0.0)), complex(expm1(-1e-20), 0.0))

@test_approx_eq expm1(complex( 1e-20, 1e-10)) complex( 5e-21, 1e-10)
@test_approx_eq expm1(complex( 1e-20,-1e-10)) complex( 5e-21,-1e-10)
@test_approx_eq expm1(complex(-1e-20, 1e-10)) complex(-1.5e-20, 1e-10)
@test_approx_eq expm1(complex(-1e-20,-1e-10)) complex(-1.5e-20,-1e-10)

@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
@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


# ^ (cpow)
# equivalent to exp(y*log(x))
# except for 0^0?
Expand Down