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: 11 on 8 virtual cores

julia> using LinearAlgebra, DynamicPolynomials

julia> function f(n)
         @PolyVar a b c d e
         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),
         )
       end
f (generic function with 1 method)

julia> const m15 = f(15);

julia> const m16 = f(16);

julia> @time det(m15);
  1.945673 seconds (45.22 M allocations: 2.261 GiB, 20.60% gc time, 4.02% compilation time)

julia> @time det(m15);
  1.991062 seconds (45.22 M allocations: 2.261 GiB, 23.74% gc time)

julia> @time det(m16);
  4.596664 seconds (106.67 M allocations: 5.324 GiB, 22.65% gc time)

julia> @time det(m16);
  4.648503 seconds (106.67 M allocations: 5.324 GiB, 22.66% gc time)
```

The above REPL session is with this commit applied, and with all other
recent PRs of mine applied, to MultivariatePolynomials.jl,
DynamicPolynomials.jl, and MutableArithmetics.jl. The same computation
with MultivariatePolynomials v0.5.3 ran for multiple minutes before I
decided to just kill it.

Depends on #285.

Fixes #281.

Author: nsajko

src/det.jl

@@ -1,13 +1,72 @@
-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}
+    get_elem = let A = A
+        (i, j) -> LinearAlgebra.det(A[i, j])
+    end
+    r = first(size(A))
+    if 1 < r
+        total = LinearAlgebra.det(zero(T))
+        for i in Base.OneTo(r)
+            a = get_elem(i, 1)
+            if !iszero(a)
+                ii = Base.OneTo(r) .!= i
+                jj = 2:r
+                B = A[ii, jj]
+                x = det_impl(B)
+                x = MA.operate!!(*, x, a)
+                if iseven(i)
+                    x = MA.operate!!(-, x)
+                end
+                total = MA.operate!!(+, total, x)
+            end
+        end
+        total
+    elseif isone(r)
+        MA.copy_if_mutable(get_elem(1, 1))
+    else
+        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
-        return M[1, 1]
+        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
+
+function promote_if_narrow(m::AbstractMatrix{<:NarrowIntegerPolynomialLike})
+    return map((p -> polynomial(p, BigInt)), m)
+end
+
+# 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))
+
+function LinearAlgebra.det(m::AbstractMatrix{<:AbstractPolynomialLike})
+    return det_impl_outer(promote_if_necessary(m))
+end

test/det.jl

@@ -1,7 +1,14 @@
 @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 in (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