Merge pull request #262 from TsurHerman/pull-request/fb369658

More  consistent treatment of unsigned types in det and inv

Author: c42f

src/det.jl

@@ -1,49 +1,26 @@
-@inline det(A::StaticMatrix) = _det(Size(A), A)
-@inline logdet(A::StaticMatrix) = _logdet(Size(A), A)
+@inline function det(A::StaticMatrix)
+    T = eltype(A)
+    S = typeof((one(T)*zero(T) + zero(T))/one(T))
+    _det(Size(A),A,S)
+end
 
-@inline _det(::Size{(1,1)}, A::StaticMatrix) = @inbounds return A[1]
+@inline logdet(A::StaticMatrix) = log(det(A))
 
-@inline function _det(::Size{(2,2)}, A::StaticMatrix)
-    @inbounds return A[1]*A[4] - A[3]*A[2]
-end
+@inline _det(::Size{(1,1)}, A::StaticMatrix,S::Type) = @inbounds return convert(S,A[1])
 
-@inline function _det(::Size{(2,2)}, A::StaticMatrix{<:Any, <:Any, <:Unsigned})
-    @inbounds return Signed(A[1]*A[4]) - Signed(A[3]*A[2])
+@inline function _det(::Size{(2,2)}, A::StaticMatrix, S::Type)
+    A = similar_type(A,S)(A)
+    @inbounds return A[1]*A[4] - A[3]*A[2]
 end
 
-@inline function _det(::Size{(3,3)}, A::StaticMatrix)
+@inline function _det(::Size{(3,3)}, A::StaticMatrix, S::Type)
+    A = similar_type(A,S)(A)
     @inbounds x0 = SVector(A[1], A[2], A[3])
     @inbounds x1 = SVector(A[4], A[5], A[6])
     @inbounds x2 = SVector(A[7], A[8], A[9])
     return vecdot(x0, cross(x1, x2))
 end
 
-@inline function _det(::Size{(3,3)}, A::StaticMatrix{<:Any, <:Any, <:Unsigned})
-    @inbounds x0 = SVector(Signed(A[1]), Signed(A[2]), Signed(A[3]))
-    @inbounds x1 = SVector(Signed(A[4]), Signed(A[5]), Signed(A[6]))
-    @inbounds x2 = SVector(Signed(A[7]), Signed(A[8]), Signed(A[9]))
-    return vecdot(x0, cross(x1, x2))
+@inline function _det(::Size, A::StaticMatrix,::Type)
+    return det(Matrix(A))
 end
-
-@inline _logdet(S::Union{Size{(1,1)},Size{(2,2)},Size{(3,3)}}, A::StaticMatrix) = log(_det(S, A))
-
-for (symb, f) in [(:_det, :det), (:_logdet, :logdet)]
-    eval(quote
-        @generated function $symb{S}(::Size{S}, A::StaticMatrix)
-            if S[1] != S[2]
-                throw(DimensionMismatch("matrix is not square"))
-            end
-            return quote # Implementation from Base
-                @_inline_meta
-                T = eltype(A)
-                T2 = typeof((one(T)*zero(T) + zero(T))/one(T))
-                if istriu(A) || istril(A)
-                    return convert(T2, $($f)(UpperTriangular(A))) # Is this a Julia bug that a convert is not type stable??
-                end
-                AA = convert(Array{T2}, A)
-                return $($f)(lufact(AA))
-            end
-        end
-    end)    
-end
-

src/inv.jl

@@ -16,9 +16,9 @@ end
     S = typeof((one(T)*zero(T) + zero(T))/one(T))
     newtype = similar_type(A, S)
 
-    @inbounds x0 = SVector(A[1], A[2], A[3])
-    @inbounds x1 = SVector(A[4], A[5], A[6])
-    @inbounds x2 = SVector(A[7], A[8], A[9])
+    @inbounds x0 = SVector{3,S}(A[1], A[2], A[3])
+    @inbounds x1 = SVector{3,S}(A[4], A[5], A[6])
+    @inbounds x2 = SVector{3,S}(A[7], A[8], A[9])
 
     y0 = cross(x1,x2)
     d  = vecdot(x0, y0)

test/det.jl

@@ -1,21 +1,24 @@
 @testset "Determinant" begin
-    @test det(@SMatrix [1]) === 1
-    @test logdet(@SMatrix [1]) === 0.0
-    @test det(@SMatrix [0 1; 1 0]) === -1
+    @test det(@SMatrix [1]) == 1
+    @test logdet(@SMatrix [1]) == 0.0
+    @test det(@SMatrix [0 1; 1 0]) == -1
     @test logdet(@SMatrix Complex{Float64}[0 1; 1 0]) == log(det(@SMatrix Complex{Float64}[0 1; 1 0]))
-    
-    @test det(@SMatrix [0 1 0; 1 0 0; 0 0 1]) === -1
+
+    @test det(@SMatrix [0 1 0; 1 0 0; 0 0 1]) == -1
     m = randn(Float64, 4,4)
     @test det(SMatrix{4,4}(m)) ≈ det(m)
     #triu/tril
-    @test det(@SMatrix [1 2; 0 3]) === 3
-    @test det(@SMatrix [1 2 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) === 400.0
+    @test det(@SMatrix [1 2; 0 3]) == 3
+    @test det(@SMatrix [1 2 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) == 400.0
     @test logdet(@SMatrix [1 2 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) ≈ log(400.0)
     @test @inferred(det(ones(SMatrix{10,10,Complex{Float64}}))) == 0
 
-    # Unsigned specializations
-    @test det(@SMatrix [0x00 0x01; 0x01 0x00])::Int8 == -1
-    @test det(@SMatrix [0x00 0x01 0x00; 0x01 0x00 0x00; 0x00 0x00 0x01])::Int8 == -1
+    # Unsigned specializations , compare to Base
+    M = @SMatrix [1 2 3 4; 200 5 6 7; 0 0 8 9; 0 0 0 10]
+    for sz in (2,3,4), typ in (UInt8,UInt16,UInt32,UInt64)
+        Mtag = SMatrix{sz,sz,typ}(M[1:sz,1:sz])
+        @test det(Mtag) == det(Array(Mtag))
+    end
 
     @test_throws DimensionMismatch det(@SMatrix [0; 1])
 end