diff --git a/src/det.jl b/src/det.jl index 66cd8dce..19736025 100644 --- a/src/det.jl +++ b/src/det.jl @@ -1,13 +1,64 @@ -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 + +# 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 + # TODO: use MutableArithmetics.jl for performance + total = (LinearAlgebra.det ∘ zero)(T) + for i ∈ Base.OneTo(r) + a = LinearAlgebra.det(A[i, 1]) + if !iszero(a) + ii = Base.OneTo(r) .!= i + jj = 2:r + B = A[ii, jj] + x = a * det_impl(B) + if iseven(i) + total -= x + else + total += x + end + end + end + total + elseif isone(r) + LinearAlgebra.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 det_impl_outer(m::AbstractMatrix{T}) where {T} + if 0 < LinearAlgebra.checksquare(m) + (det_impl ∘ collect_if_not_already_matrix)(m) + else + (LinearAlgebra.det ∘ one)(T) + end +end + +# Determinants of narrow integer type: `LinearAlgebra` seems to +# promote these to `Float64` to prevent them from overflowing. We +# instead promote to `BigInt` to keep things exact. In the case of +# `Bool` we also need to promote for type stability. + +const DeterminantNarrowIntegerTypes = Union{ + Bool, UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64, + UInt128, Int128, +} + +const NarrowIntegerPolynomialLike = + AbstractPolynomialLike{T} where {T<:DeterminantNarrowIntegerTypes} + +promote_if_narrow(m::AbstractMatrix{<:AbstractPolynomialLike}) = m + +# TODO: avoid the unnecessary multiplication +promote_if_narrow(m::AbstractMatrix{<:NarrowIntegerPolynomialLike}) = + m * one(BigInt) + +LinearAlgebra.det(m::AbstractMatrix{<:AbstractPolynomialLike}) = + (det_impl_outer ∘ promote_if_narrow)(m)