Fix matrix trig functions for matrices of integers (#24180)

* Convert integers to float in matrix trig functions * Add tests * Dots!
JuliaLang · Oct 23, 2017 · 8798936 · 8798936
1 parent 9874daa
commit 8798936
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 12 deletions.
diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
@@ -720,14 +720,16 @@ function cos(A::AbstractMatrix{<:Real})
     if issymmetric(A)
         return copytri!(parent(cos(Symmetric(A))), 'U')
     end
-    return real(exp!(im*A))
+    T = complex(float(eltype(A)))
+    return real(exp!(T.(im .* A)))
 end
 function cos(A::AbstractMatrix{<:Complex})
     if ishermitian(A)
         return copytri!(parent(cos(Hermitian(A))), 'U', true)
     end
-    X = exp!(im*A)
-    X .= (X .+ exp!(-im*A)) ./ 2
+    T = complex(float(eltype(A)))
+    X = exp!(T.(im .* A))
+    @. X = (X + $exp!(T(-im*A))) / 2
     return X
 end
 
@@ -751,14 +753,16 @@ function sin(A::AbstractMatrix{<:Real})
     if issymmetric(A)
         return copytri!(parent(sin(Symmetric(A))), 'U')
     end
-    return imag(exp!(im*A))
+    T = complex(float(eltype(A)))
+    return imag(exp!(T.(im .* A)))
 end
 function sin(A::AbstractMatrix{<:Complex})
     if ishermitian(A)
         return copytri!(parent(sin(Hermitian(A))), 'U', true)
     end
-    X = exp!(im*A)
-    Y = exp!(-im*A)
+    T = complex(float(eltype(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 +797,8 @@ function sincos(A::AbstractMatrix{<:Real})
         cosA = copytri!(parent(symcosA), 'U')
         return sinA, cosA
     end
-    c, s = reim(exp!(im*A))
+    T = complex(float(eltype(A)))
+    c, s = reim(exp!(T.(im .* A)))
     return s, c
 end
 function sincos(A::AbstractMatrix{<:Complex})
@@ -803,8 +808,9 @@ function sincos(A::AbstractMatrix{<:Complex})
         cosA = copytri!(parent(hermcosA), 'U', true)
         return sinA, cosA
     end
-    X = exp!(im*A)
-    Y = exp!(-im*A)
+    T = complex(float(eltype(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))
@@ -848,7 +854,7 @@ function cosh(A::AbstractMatrix)
         return copytri!(parent(cosh(Hermitian(A))), 'U', true)
     end
     X = exp(A)
-    X .= (X .+ exp!(-A)) ./ 2
+    @. X = (X + $exp!(float(-A))) / 2
     return X
 end
 
@@ -862,7 +868,7 @@ function sinh(A::AbstractMatrix)
         return copytri!(parent(sinh(Hermitian(A))), 'U', true)
     end
     X = exp(A)
-    X .= (X .- exp!(-A)) ./ 2
+    @. X = (X - $exp!(float(-A))) / 2
     return X
 end
 
@@ -876,7 +882,7 @@ function tanh(A::AbstractMatrix)
         return copytri!(parent(tanh(Hermitian(A))), 'U', true)
     end
     X = exp(A)
-    Y = exp!(-A)
+    Y = exp!(float.(.-A))
     @inbounds for i in eachindex(X)
         x, y = X[i], Y[i]
         X[i] = x - y

diff --git 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])