Improve performance of derangement/subfactorial with iterative implementation (#146)

* Improve performance of `derangement`/`subfactorial` with iterative implementation

Use the recursive formula
  !n = (n-1) * (!(n-1) + !(n-2))
presented here: https://en.wikipedia.org/wiki/Derangement#Counting_derangements

* Improve performance of `derangement`/`subfactorial` with simpler recursive formula and inplace computations

Use the simpler formula

!n = n * !(n-1) + (-1)^n

and use inplace operations on `BigInt`s to avoid allocations.

* Update test/factorials.jl

---------

Co-authored-by: Chengyu Han <[email protected]>

Author: FedericoStra

src/factorials.jl

@@ -38,9 +38,23 @@ Base.factorial(n::Integer, k::Integer) = factorial(promote(n, k)...)
 Compute the number of permutations of `n` with no fixed points, also known as the
 subfactorial. An alias `subfactorial` for this function is provided for convenience.
 """
-function derangement(sn::Integer)
-    n = BigInt(sn)
-    return numerator(factorial(n) * sum([(-1)^k // factorial(k) for k = 0:n]))
+function derangement(n::Integer)
+    if n < 0
+        throw(DomainError(n, "n must be nonnegative"))
+    elseif n <= 1
+        return BigInt(1-n)
+    end
+    d = BigInt(0)
+    for i in 2:n
+        # d = i * d + (iseven(i) ? 1 : -1)
+        Base.GMP.MPZ.mul_ui!(d, i)
+        if iseven(i)
+            Base.GMP.MPZ.add_ui!(d, 1)
+        else
+            Base.GMP.MPZ.sub_ui!(d, 1)
+        end
+    end
+    return d
 end
 const subfactorial = derangement
 

test/factorials.jl

@@ -10,7 +10,10 @@
 
     # derangement
     @test derangement(4) == subfactorial(4) == 9
+    @test derangement(0) == 1
+    @test derangement(1) == 0
     @test derangement(24) == parse(BigInt, "228250211305338670494289")
+    @test_throws DomainError derangement(-1)
 
     # partialderangement
     @test partialderangement(7, 3) == 315