New permanent algorithm

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearAlgebraX"
 uuid = "9b3f67b0-2d00-526e-9884-9e4938f8fb88"
-version = "0.2.5"
+version = "0.2.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/perm.jl

@@ -1,16 +1,43 @@
 export permanent
+
+_enlarge(x::T) where {T<:Union{IntegerX,RationalX}} = big(x)
+_enlarge(x) = x
+
 """
     permanent(A)
 
-Compute the permanent of the square matrix `A`. Works best if `A`
-is relatively sparse. For a dense matrix, this will be slow.
+Compute the permanent of the square matrix `A`.
 """
-function permanent(A::AbstractMatrix{T}) where {T<:Union{IntegerX,RationalX}}
+function permanent(A::AbstractMatrix{T}) where {T}     # code by Daniel Scheinerman
+    A = _enlarge.(A)
+    m, n = size(A)
+    s = _enlarge(zero(T)) #to accumulate the permanent
+    x = UInt64(0) #gray code state
+    v = zeros(Int64, n) #current subset sum
+    for i = 1:2^n-1
+        t = trailing_zeros(i) #next bit flip
+        b = (x >> t) & 1 #if 1 then add row, else subtract
+        x = xor(x, 1 << t)
+        if b == 0
+            v += A[t+1, :]
+        else
+            v -= A[t+1, :]
+        end
+        s += (-1)^(i & 1) * prod(v)
+    end
+    s *= (-1)^n
+    return s
+end
+
+
+# Old code available as LinearAlgebraX.old_permanent
+
+function old_permanent(A::AbstractMatrix{T}) where {T<:Union{IntegerX,RationalX}}
     A = big.(A)
     return perm_work(A)
 end
 
-function permanent(A::AbstractMatrix{T})::T where {T}
+function old_permanent(A::AbstractMatrix{T})::T where {T}
     perm_work(A)
 end