From e5056f53967c11f9e8fb9b05a03b22d18b0bc400 Mon Sep 17 00:00:00 2001 From: Evey Dee Date: Tue, 17 Oct 2017 15:37:33 +0200 Subject: [PATCH 1/3] Convert integers to float in matrix trig functions --- base/linalg/dense.jl | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index f0995853d3b6b..065e81d7a413c 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -718,14 +718,16 @@ function cos(A::AbstractMatrix{<:Real}) if issymmetric(A) return full(cos(Symmetric(A))) end - return real(exp!(im*A)) + T = complex(float(eltype(A))) + return real(exp!(convert(AbstractArray{T}, im*A))) end function cos(A::AbstractMatrix{<:Complex}) if ishermitian(A) return full(cos(Hermitian(A))) end - X = exp!(im*A) - X .= (X .+ exp!(-im*A)) ./ 2 + T = complex(float(eltype(A))) + X = exp!(convert(AbstractArray{T}, im*A)) + X .= (X .+ exp!(convert(AbstractArray{T}, -im*A))) ./ 2 return X end @@ -749,14 +751,16 @@ function sin(A::AbstractMatrix{<:Real}) if issymmetric(A) return full(sin(Symmetric(A))) end - return imag(exp!(im*A)) + T = complex(float(eltype(A))) + return imag(exp!(convert(AbstractArray{T}, im*A))) end function sin(A::AbstractMatrix{<:Complex}) if ishermitian(A) return full(sin(Hermitian(A))) end - X = exp!(im*A) - Y = exp!(-im*A) + T = complex(float(eltype(A))) + X = exp!(convert(AbstractArray{T}, im*A)) + Y = exp!(convert(AbstractArray{T}, -im*A)) @inbounds for i in eachindex(X) x, y = X[i]/2, Y[i]/2 X[i] = Complex(imag(x)-imag(y), real(y)-real(x)) @@ -788,15 +792,17 @@ function sincos(A::AbstractMatrix{<:Real}) if issymmetric(A) return full.(sincos(Symmetric(A))) end - c, s = reim(exp!(im*A)) + T = complex(float(eltype(A))) + c, s = reim(exp!(convert(AbstractArray{T}, im*A))) return s, c end function sincos(A::AbstractMatrix{<:Complex}) if ishermitian(A) return full.(sincos(Hermitian(A))) end - X = exp!(im*A) - Y = exp!(-im*A) + T = complex(float(eltype(A))) + X = exp!(convert(AbstractArray{T}, im*A)) + Y = exp!(convert(AbstractArray{T}, -im*A)) @inbounds for i in eachindex(X) x, y = X[i]/2, Y[i]/2 X[i] = Complex(imag(x)-imag(y), real(y)-real(x)) @@ -839,8 +845,9 @@ function cosh(A::AbstractMatrix) if ishermitian(A) return full(cosh(Hermitian(A))) end + T = float(eltype(A)) X = exp(A) - X .= (X .+ exp!(-A)) ./ 2 + X .= (X .+ exp!(convert(AbstractArray{T}, -A))) ./ 2 return X end @@ -853,8 +860,9 @@ function sinh(A::AbstractMatrix) if ishermitian(A) return full(sinh(Hermitian(A))) end + T = float(eltype(A)) X = exp(A) - X .= (X .- exp!(-A)) ./ 2 + X .= (X .- exp!(convert(AbstractArray{T}, -A))) ./ 2 return X end @@ -867,8 +875,9 @@ function tanh(A::AbstractMatrix) if ishermitian(A) return full(tanh(Hermitian(A))) end + T = float(eltype(A)) X = exp(A) - Y = exp!(-A) + Y = exp!(convert(AbstractArray{T}, -A)) @inbounds for i in eachindex(X) x, y = X[i], Y[i] X[i] = x - y From 8008ec6fb53bf2a9dc08a7e45f9b7eaee9ee728a Mon Sep 17 00:00:00 2001 From: Evey Dee Date: Tue, 17 Oct 2017 15:55:26 +0200 Subject: [PATCH 2/3] Add tests --- test/linalg/dense.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index 1ec876a602bb9..69b38d09ec892 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -508,6 +508,27 @@ end @test sinh(A5) ≈ 0.5 * (exp(A5) - exp(-A5)) end + @testset "Additional tests for $elty" for elty in (Int32, Int64, Complex{Int32}, Complex{Int64}) + A1 = convert(Matrix{elty}, [1 2; 3 4]) + A2 = convert(Matrix{elty}, [1 2; 2 1]) + + cosA1 = convert(Matrix{float(elty)}, [0.855423165077998 -0.11087638101074865; + -0.16631457151612294 0.689108593561875]) + cosA2 = convert(Matrix{float(elty)}, [-0.22484509536615283 -0.7651474012342925; + -0.7651474012342925 -0.22484509536615283]) + + @test cos(A1) ≈ cosA1 + @test cos(A2) ≈ cosA2 + + sinA1 = convert(Matrix{float(elty)}, [-0.46558148631373036 -0.14842445991317652; + -0.22263668986976476 -0.6882181761834951]) + sinA2 = convert(Matrix{float(elty)}, [-0.3501754883740146 0.4912954964338818; + 0.4912954964338818 -0.3501754883740146]) + + @test sin(A1) ≈ sinA1 + @test sin(A2) ≈ sinA2 + end + @testset "Inverse functions for $elty" for elty in (Float32, Float64) A1 = convert(Matrix{elty}, [0.244637 -0.63578; 0.22002 0.189026]) From 1ccecc9d17a7767508ebbbda7caaa001c1a67db0 Mon Sep 17 00:00:00 2001 From: Evey Dee Date: Mon, 23 Oct 2017 17:06:48 +0200 Subject: [PATCH 3/3] Dots! --- base/linalg/dense.jl | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index 065e81d7a413c..16edf625ce8b7 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -719,15 +719,15 @@ function cos(A::AbstractMatrix{<:Real}) return full(cos(Symmetric(A))) end T = complex(float(eltype(A))) - return real(exp!(convert(AbstractArray{T}, im*A))) + return real(exp!(T.(im .* A))) end function cos(A::AbstractMatrix{<:Complex}) if ishermitian(A) return full(cos(Hermitian(A))) end T = complex(float(eltype(A))) - X = exp!(convert(AbstractArray{T}, im*A)) - X .= (X .+ exp!(convert(AbstractArray{T}, -im*A))) ./ 2 + X = exp!(T.(im .* A)) + @. X = (X + $exp!(T(-im*A))) / 2 return X end @@ -752,15 +752,15 @@ function sin(A::AbstractMatrix{<:Real}) return full(sin(Symmetric(A))) end T = complex(float(eltype(A))) - return imag(exp!(convert(AbstractArray{T}, im*A))) + return imag(exp!(T.(im .* A))) end function sin(A::AbstractMatrix{<:Complex}) if ishermitian(A) return full(sin(Hermitian(A))) end T = complex(float(eltype(A))) - X = exp!(convert(AbstractArray{T}, im*A)) - Y = exp!(convert(AbstractArray{T}, -im*A)) + X = exp!(T.(im .* A)) + Y = exp!(T.(.-im .* A)) @inbounds for i in eachindex(X) x, y = X[i]/2, Y[i]/2 X[i] = Complex(imag(x)-imag(y), real(y)-real(x)) @@ -793,7 +793,7 @@ function sincos(A::AbstractMatrix{<:Real}) return full.(sincos(Symmetric(A))) end T = complex(float(eltype(A))) - c, s = reim(exp!(convert(AbstractArray{T}, im*A))) + c, s = reim(exp!(T.(im .* A))) return s, c end function sincos(A::AbstractMatrix{<:Complex}) @@ -801,8 +801,8 @@ function sincos(A::AbstractMatrix{<:Complex}) return full.(sincos(Hermitian(A))) end T = complex(float(eltype(A))) - X = exp!(convert(AbstractArray{T}, im*A)) - Y = exp!(convert(AbstractArray{T}, -im*A)) + X = exp!(T.(im .* A)) + Y = exp!(T.(.-im .* A)) @inbounds for i in eachindex(X) x, y = X[i]/2, Y[i]/2 X[i] = Complex(imag(x)-imag(y), real(y)-real(x)) @@ -845,9 +845,8 @@ function cosh(A::AbstractMatrix) if ishermitian(A) return full(cosh(Hermitian(A))) end - T = float(eltype(A)) X = exp(A) - X .= (X .+ exp!(convert(AbstractArray{T}, -A))) ./ 2 + @. X = (X + $exp!(float(-A))) / 2 return X end @@ -860,9 +859,8 @@ function sinh(A::AbstractMatrix) if ishermitian(A) return full(sinh(Hermitian(A))) end - T = float(eltype(A)) X = exp(A) - X .= (X .- exp!(convert(AbstractArray{T}, -A))) ./ 2 + @. X = (X - $exp!(float(-A))) / 2 return X end @@ -875,9 +873,8 @@ function tanh(A::AbstractMatrix) if ishermitian(A) return full(tanh(Hermitian(A))) end - T = float(eltype(A)) X = exp(A) - Y = exp!(convert(AbstractArray{T}, -A)) + Y = exp!(float.(.-A)) @inbounds for i in eachindex(X) x, y = X[i], Y[i] X[i] = x - y