diff --git a/src/det.jl b/src/det.jl index 66cd8dce..c2593b22 100644 --- a/src/det.jl +++ b/src/det.jl @@ -1,13 +1,42 @@ -function LinearAlgebra.det(M::Matrix{<:AbstractPolynomialLike}) - m = size(M)[1] - if m > 2 - return sum( - (-1)^(i - 1) * M[i, 1] * LinearAlgebra.det(M[1:end.!=i, 2:end]) for - i in 1:m - ) - elseif m == 2 - return M[1, 1] * M[2, 2] - M[2, 1] * M[1, 2] +# Scalar determinant, used for recursive computation of the determinant +LinearAlgebra.det(p::AbstractPolynomialLike{<:Number}) = p + +element_det(e) = LinearAlgebra.det(e)*one(Int8) + +# Matrix determinant by cofactor expansion, adapted from +# `LinearAlgebraX.cofactor_det`. + +function det_impl(A::AbstractMatrix{T}) where {T} + r = (first ∘ size)(A) + if 1 < r + total = (element_det ∘ zero)(T) + for i ∈ Base.OneTo(r) + a = element_det(A[i, 1]) + if !iszero(a) + ii = Base.OneTo(r) .!= i + jj = 2:r + B = A[ii, jj] + sign = let o = one(Int8) + iseven(i) ? -o : o + end + total += sign * a * det_impl(B) + end + end + total + elseif isone(r) + element_det(A[1, 1]) else - return M[1, 1] + error("unexpected") + end +end + +collect_if_not_already_matrix(m::Matrix) = m +collect_if_not_already_matrix(m::AbstractMatrix) = collect(m) + +function LinearAlgebra.det(m::AbstractMatrix{<:AbstractPolynomialLike}) + if 0 < LinearAlgebra.checksquare(m) + (det_impl ∘ collect_if_not_already_matrix)(m) + else + (element_det ∘ one ∘ eltype)(m) end end