LinearAlgebra.det improvements

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 -t 8
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.11.0-DEV.972 (2023-11-23)
 _/ |\__'_|_|_|\__'_|  |  Commit 9884e447e79 (1 day old master)
|__/                   |

julia> using LinearAlgebra, DynamicPolynomials

julia> @PolyVar a b c d e
(a, b, c, d, e)

julia> const m = diagm(
         -2 => fill(a, 14), -1 => fill(b, 15),
         0  => fill(c, 16),
         2  => fill(e, 14), 1  => fill(d, 15))
16×16 Matrix{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Int64}}:
 c  d  e  0  0  0  0  0  0  0  0  0  0  0  0  0
 b  c  d  e  0  0  0  0  0  0  0  0  0  0  0  0
 a  b  c  d  e  0  0  0  0  0  0  0  0  0  0  0
 0  a  b  c  d  e  0  0  0  0  0  0  0  0  0  0
 0  0  a  b  c  d  e  0  0  0  0  0  0  0  0  0
 0  0  0  a  b  c  d  e  0  0  0  0  0  0  0  0
 0  0  0  0  a  b  c  d  e  0  0  0  0  0  0  0
 0  0  0  0  0  a  b  c  d  e  0  0  0  0  0  0
 0  0  0  0  0  0  a  b  c  d  e  0  0  0  0  0
 0  0  0  0  0  0  0  a  b  c  d  e  0  0  0  0
 0  0  0  0  0  0  0  0  a  b  c  d  e  0  0  0
 0  0  0  0  0  0  0  0  0  a  b  c  d  e  0  0
 0  0  0  0  0  0  0  0  0  0  a  b  c  d  e  0
 0  0  0  0  0  0  0  0  0  0  0  a  b  c  d  e
 0  0  0  0  0  0  0  0  0  0  0  0  a  b  c  d
 0  0  0  0  0  0  0  0  0  0  0  0  0  a  b  c

julia> @time det(m);
  9.301086 seconds (162.91 M allocations: 5.968 GiB, 9.55% gc time, 1.93% compilation time)

julia> @time det(m);
  9.197843 seconds (162.88 M allocations: 5.967 GiB, 10.97% gc time)

julia> @time det(m);
  9.628294 seconds (162.88 M allocations: 5.967 GiB, 10.74% 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,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