From 61d6384636a4fdfdd56d41ead2c3c7bd5e970702 Mon Sep 17 00:00:00 2001 From: ExpandingMan <savastio@protonmail.com> Date: Tue, 23 Oct 2018 13:10:29 -0400 Subject: [PATCH 1/6] added missing matrix exponentiation methods --- base/mathconstants.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/base/mathconstants.jl b/base/mathconstants.jl index a3d1be99becbb..721a429f1cbda 100644 --- a/base/mathconstants.jl +++ b/base/mathconstants.jl @@ -89,6 +89,7 @@ catalan # loop over types to prevent ambiguities for ^(::Number, x) for T in (AbstractIrrational, Rational, Integer, Number, Complex) Base.:^(::Irrational{:ℯ}, x::T) = exp(x) + Base.:^(::Irrational{:ℯ}, x::AbstractMatrix{T}) = exp(x) end Base.literal_pow(::typeof(^), ::Irrational{:ℯ}, ::Val{p}) where {p} = exp(p) From 9881db3885fc3ddfc1003674f7f45d0312da4735 Mon Sep 17 00:00:00 2001 From: ExpandingMan <savastio@protonmail.com> Date: Tue, 23 Oct 2018 13:22:36 -0400 Subject: [PATCH 2/6] added some tests --- test/math.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/math.jl b/test/math.jl index 85b81b6f17e84..eba799d03876c 100644 --- a/test/math.jl +++ b/test/math.jl @@ -40,6 +40,11 @@ end @test ℯ^2 == exp(2) @test ℯ^2.4 == exp(2.4) @test ℯ^(2//3) == exp(2//3) + @test ℯ^(fill(2, (4,4))) == exp(fill(2, (4,4))) + @test ℯ^(fill(2.4, (4,4))) == exp(fill(2.4, (4,4))) + # currently missing `exp` methods for matrices of Rational and Irrational + # @test ℯ^(fill(π, (4,4))) == exp(fill(π, (4,4))) + # @test ℯ^(fill(2//3, (4,4))) == exp(fill(2//3, (4,4))) @test Float16(3.0) < pi @test pi < Float16(4.0) From 0c483a8e0cf6e208b74bde48f6cc866239a4d46e Mon Sep 17 00:00:00 2001 From: ExpandingMan <savastio@protonmail.com> Date: Tue, 23 Oct 2018 15:29:59 -0400 Subject: [PATCH 3/6] moved new methods to LinearAlgebra --- base/mathconstants.jl | 1 - stdlib/LinearAlgebra/src/LinearAlgebra.jl | 3 +++ stdlib/LinearAlgebra/test/dense.jl | 4 ++++ test/math.jl | 5 ----- 4 files changed, 7 insertions(+), 6 deletions(-) diff --git a/base/mathconstants.jl b/base/mathconstants.jl index 721a429f1cbda..a3d1be99becbb 100644 --- a/base/mathconstants.jl +++ b/base/mathconstants.jl @@ -89,7 +89,6 @@ catalan # loop over types to prevent ambiguities for ^(::Number, x) for T in (AbstractIrrational, Rational, Integer, Number, Complex) Base.:^(::Irrational{:ℯ}, x::T) = exp(x) - Base.:^(::Irrational{:ℯ}, x::AbstractMatrix{T}) = exp(x) end Base.literal_pow(::typeof(^), ::Irrational{:ℯ}, ::Val{p}) where {p} = exp(p) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index c58f98dd83f31..029a649ac958c 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -340,6 +340,9 @@ rdiv!(A, B) copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A) copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A) +# matrix exponentials for ℯ +Base.:^(::Irrational{:ℯ}, A::AbstractMatrix) = exp(A) + include("adjtrans.jl") include("transpose.jl") diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl index ae935b6337b59..b9ca0bc7f4e08 100644 --- a/stdlib/LinearAlgebra/test/dense.jl +++ b/stdlib/LinearAlgebra/test/dense.jl @@ -401,6 +401,10 @@ end [4.000000000000000 -1.414213562373094 -1.414213562373095 -1.414213562373095 4.999999999999996 -0.000000000000000 0 -0.000000000000002 3.000000000000000]) + + # alias using ℯ, should be exact + @test ℯ^(fill(2, (4,4))) == exp(fill(2, (4,4))) + @test ℯ^(fill(2.4, (4,4))) == exp(fill(2.4, (4,4))) end @testset "Additional tests for $elty" for elty in (Float64, Complex{Float64}) diff --git a/test/math.jl b/test/math.jl index eba799d03876c..85b81b6f17e84 100644 --- a/test/math.jl +++ b/test/math.jl @@ -40,11 +40,6 @@ end @test ℯ^2 == exp(2) @test ℯ^2.4 == exp(2.4) @test ℯ^(2//3) == exp(2//3) - @test ℯ^(fill(2, (4,4))) == exp(fill(2, (4,4))) - @test ℯ^(fill(2.4, (4,4))) == exp(fill(2.4, (4,4))) - # currently missing `exp` methods for matrices of Rational and Irrational - # @test ℯ^(fill(π, (4,4))) == exp(fill(π, (4,4))) - # @test ℯ^(fill(2//3, (4,4))) == exp(fill(2//3, (4,4))) @test Float16(3.0) < pi @test pi < Float16(4.0) From 72266ad9bbb306b7be1d5deaf43b8c95468bc7e3 Mon Sep 17 00:00:00 2001 From: ExpandingMan <savastio@protonmail.com> Date: Wed, 24 Oct 2018 09:58:15 -0400 Subject: [PATCH 4/6] added general method and cleaned up tests --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 4 +++- stdlib/LinearAlgebra/test/dense.jl | 11 +++++++---- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 029a649ac958c..4e939ee339af8 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -340,7 +340,9 @@ rdiv!(A, B) copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A) copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A) -# matrix exponentials for ℯ +# matrix exponentials +Base.:^(b::Number, A::AbstractMatrix) = exp(log(b)*A) +# method for ℯ to explicitly elide the log(b) multiplication Base.:^(::Irrational{:ℯ}, A::AbstractMatrix) = exp(A) include("adjtrans.jl") diff --git a/stdlib/LinearAlgebra/test/dense.jl b/stdlib/LinearAlgebra/test/dense.jl index b9ca0bc7f4e08..32b217b6fdc53 100644 --- a/stdlib/LinearAlgebra/test/dense.jl +++ b/stdlib/LinearAlgebra/test/dense.jl @@ -401,10 +401,6 @@ end [4.000000000000000 -1.414213562373094 -1.414213562373095 -1.414213562373095 4.999999999999996 -0.000000000000000 0 -0.000000000000002 3.000000000000000]) - - # alias using ℯ, should be exact - @test ℯ^(fill(2, (4,4))) == exp(fill(2, (4,4))) - @test ℯ^(fill(2.4, (4,4))) == exp(fill(2.4, (4,4))) end @testset "Additional tests for $elty" for elty in (Float64, Complex{Float64}) @@ -432,6 +428,13 @@ end end end + @testset "^ tests" for elty in (Float32, Float64, ComplexF32, ComplexF64, Int32, Int64) + # should all be exact as the lhs functions are simple aliases + @test ℯ^(fill(elty(2), (4,4))) == exp(fill(elty(2), (4,4))) + @test 2^(fill(elty(2), (4,4))) == exp(log(2)*fill(elty(2), (4,4))) + @test 2.0^(fill(elty(2), (4,4))) == exp(log(2.0)*fill(elty(2), (4,4))) + end + A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1] @test exp(log(A8)) ≈ A8 end From 9def8c73a7bc2e6ea9e6943c86d683f1e6ad7ce8 Mon Sep 17 00:00:00 2001 From: ExpandingMan <savastio@protonmail.com> Date: Wed, 24 Oct 2018 10:30:36 -0400 Subject: [PATCH 5/6] moved methods --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 5 ----- stdlib/LinearAlgebra/src/dense.jl | 4 ++++ 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 4e939ee339af8..c58f98dd83f31 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -340,11 +340,6 @@ rdiv!(A, B) copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A) copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A) -# matrix exponentials -Base.:^(b::Number, A::AbstractMatrix) = exp(log(b)*A) -# method for ℯ to explicitly elide the log(b) multiplication -Base.:^(::Irrational{:ℯ}, A::AbstractMatrix) = exp(A) - include("adjtrans.jl") include("transpose.jl") diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 81a50391b91ce..d3fd085e6a473 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -489,6 +489,10 @@ julia> exp(A) exp(A::StridedMatrix{<:BlasFloat}) = exp!(copy(A)) exp(A::StridedMatrix{<:Union{Integer,Complex{<:Integer}}}) = exp!(float.(A)) +Base.:^(b::Number, A::AbstractMatrix) = exp(log(b)*A) +# method for ℯ to explicitly elide the log(b) multiplication +Base.:^(::Irrational{:ℯ}, A::AbstractMatrix) = exp(A) + ## Destructive matrix exponential using algorithm from Higham, 2008, ## "Functions of Matrices: Theory and Computation", SIAM function exp!(A::StridedMatrix{T}) where T<:BlasFloat From 6b9f486161ca0ddf6caed5d735d81a9bd8017929 Mon Sep 17 00:00:00 2001 From: ExpandingMan <ExpandingMan@users.noreply.github.com> Date: Mon, 29 Oct 2018 18:49:07 -0400 Subject: [PATCH 6/6] eliminate spurious allocation --- stdlib/LinearAlgebra/src/dense.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index d3fd085e6a473..e83164c01216a 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -489,7 +489,7 @@ julia> exp(A) exp(A::StridedMatrix{<:BlasFloat}) = exp!(copy(A)) exp(A::StridedMatrix{<:Union{Integer,Complex{<:Integer}}}) = exp!(float.(A)) -Base.:^(b::Number, A::AbstractMatrix) = exp(log(b)*A) +Base.:^(b::Number, A::AbstractMatrix) = exp!(log(b)*A) # method for ℯ to explicitly elide the log(b) multiplication Base.:^(::Irrational{:ℯ}, A::AbstractMatrix) = exp(A)