LinearAlgebra.det improvements

Performance improvements, support for more types.

Still broken for `LinearAlgebra.Symmetric` polynomial matrices,
producing a `MethodError` because of a missing `oneunit` method. This,
however, seems like a separate matter that would better be addressed
by a separate pull request.

Performance comparison:

```julia-repl
julia> versioninfo()
Julia Version 1.11.0-DEV.972
Commit 9884e447e79 (2023-11-23 16:16 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 1 on 8 virtual cores

julia> using LinearAlgebra, DynamicPolynomials

julia> @PolyVar a b c d e;

julia> f(n) = diagm(
          -2 => fill(a, n - 2), -1 => fill(b, n - 1), 0 => fill(c, n),
           2 => fill(e, n - 2),  1 => fill(d, n - 1),
       )
f (generic function with 1 method)

julia> const m15 = f(15);

julia> const m16 = f(16);

julia> @time det(m15);
  2.834004 seconds (47.98 M allocations: 1.832 GiB, 16.66% gc time, 10.61% compilation time)

julia> @time det(m15);
  2.647003 seconds (47.94 M allocations: 1.830 GiB, 18.63% gc time)

julia> @time det(m16);
  6.441887 seconds (112.61 M allocations: 4.299 GiB, 20.25% gc time)

julia> @time det(m16);
  6.390883 seconds (112.61 M allocations: 4.299 GiB, 19.33% gc time)
```

The above REPL session is with this commit applied. The same
computation with MultivariatePolynomials v0.5.3 ran for multiple
minutes before I decided to just kill it.

Fixes #281

Author: nsajko

src/det.jl

@@ -1,13 +1,70 @@
-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 NarrowIntegerTypes = Union{
+    Bool, UInt8, Int8, UInt16, Int16, UInt32, Int32, UInt64, Int64,
+    UInt128, Int128,
+}
+
+const NarrowIntegerPolynomialLike =
+  AbstractPolynomialLike{T} where {T<:NarrowIntegerTypes}
+
+promote_if_narrow(m::AbstractMatrix{<:AbstractPolynomialLike}) = m
+
+promote_if_narrow(m::AbstractMatrix{<:NarrowIntegerPolynomialLike}) =
+  map((p -> polynomial(p, BigInt)), m)
+
+# For type stability, we want to promote termlikes to polynomiallikes
+# before attempting to calculate the determinant.
+promote_if_termlike(m::AbstractMatrix{<:AbstractPolynomialLike}) = m
+promote_if_termlike(m::AbstractMatrix{<:AbstractTermLike}) = map(polynomial, m)
+
+promote_if_necessary(m) = promote_if_termlike(promote_if_narrow(m))
+
+LinearAlgebra.det(m::AbstractMatrix{<:AbstractPolynomialLike}) =
+    det_impl_outer(promote_if_necessary(m))

test/det.jl

@@ -1,7 +1,13 @@
 @testset "Det" begin
     Mod.@polyvar x y
 
-    @test det([x 1; y 1]) == x - y
-    @test det([x 1; x 1]) == 0
-    @test det([x 1 1 1; x y 1 2; 0 0 0 1; x 0 y 0]) == -x * y^2 + 2 * x * y - x
+    @testset "zero-one $T" for T ∈ (Bool, Int, Float64, BigInt, BigFloat)
+        z = zero(T)
+        o = one(T)
+        @test (@inferred det([x o; y o])) == x - y
+        @test (@inferred det([x o; x o])) == 0
+        @test (@inferred det([x+y y; z y])) == (x+y)*y
+    end
+
+    @test (@inferred det([x 1 1 1; x y 1 2; 0 0 0 1; x 0 y 0])) == -x * y^2 + 2 * x * y - x
 end