From ddd501b092abc86d9c39701055b89020ac5b1043 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Fri, 24 Nov 2023 12:41:51 +0100
Subject: [PATCH] `LinearAlgebra.det` improvements
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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> 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> @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.849449 seconds (47.98 M allocations: 1.832 GiB, 16.50% gc time, 11.14% compilation time)

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

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

julia> @time det(m16);
  6.737904 seconds (112.61 M allocations: 4.299 GiB, 22.30% 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
---
 src/det.jl | 71 ++++++++++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 61 insertions(+), 10 deletions(-)

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)