From 1d31b075fdc815baf6e507eac1d33d2beae7f5be 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

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.
---
 src/det.jl  | 79 ++++++++++++++++++++++++++++++++++++++++++++++-------
 test/det.jl | 13 +++++++--
 2 files changed, 79 insertions(+), 13 deletions(-)

diff --git a/src/det.jl b/src/det.jl
index 66cd8dce..dd352221 100644
--- a/src/det.jl
+++ b/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
diff --git a/test/det.jl b/test/det.jl
index dad8f2fc..1ef44d48 100644
--- a/test/det.jl
+++ b/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