faster permutations, based on what rdeits suggested

Project.toml

@@ -7,7 +7,8 @@ julia = "1"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["LinearAlgebra", "Test"]
+test = ["LinearAlgebra", "OffsetArrays", "Test"]

src/permutations.jl

@@ -10,66 +10,122 @@ export
 
 
 struct Permutations{T}
-    a::T
-    t::Int
+    data::T
+    length::Int
 end
 
-Base.eltype(::Type{Permutations{T}}) where {T} = Vector{eltype(T)}
+function has_repeats(state::Vector{Int})
+    # This can be safely marked inbounds because of the type restriction in the signature.
+    # If the type restriction is ever loosened, please check safety of the `@inbounds`
+    @inbounds for outer in eachindex(state)
+        for inner in (outer+1):lastindex(state)
+            if state[outer] == state[inner]
+                return true
+            end
+        end
+    end
+    return false
+end
+
+function increment!(state::Vector{Int}, min::Int, max::Int)
+    state[end] += 1
+    for i in reverse(eachindex(state))[firstindex(state):end-1]
+        if state[i] > max
+            state[i] = min
+            state[i-1] += 1
+        end
+    end
+end
+
+function next_permutation!(state::Vector{Int}, min::Int, max::Int)
+    while true
+        increment!(state, min, max)
+        has_repeats(state) || break
+    end
+end
+
+function Base.iterate(p::Permutations, state::Vector{Int}=fill(firstindex(p.data), p.length))
+    next_permutation!(state, firstindex(p.data), lastindex(p.data))
+    if first(state) > lastindex(p.data)
+        return nothing
+    end
+    [p.data[i] for i in state], state
+end
+
+function Base.length(p::Permutations)
+    length(p.data) < p.length && return 0
+    return Int(prod(length(p.data) - p.length + 1:length(p.data)))
+end
+
+Base.eltype(p::Permutations) = Vector{eltype(p.data)}
+
+Base.IteratorSize(p::Permutations) = Base.HasLength()
 
-Base.length(p::Permutations) = (0 <= p.t <= length(p.a)) ? factorial(length(p.a), length(p.a)-p.t) : 0
 
 """
     permutations(a)
 
 Generate all permutations of an indexable object `a` in lexicographic order. Because the number of permutations
 can be very large, this function returns an iterator object.
 Use `collect(permutations(a))` to get an array of all permutations.
+Only works for `a` with defined length.
 """
-permutations(a) = Permutations(a, length(a))
+permutations(a) = permutations(a, length(a))
 
 """
     permutations(a, t)
 
 Generate all size `t` permutations of an indexable object `a`.
+Only works for `a` with defined length. 
+If `(t <= 0) || (t > length(a))`, then returns an empty vector of eltype of `a`
 """
 function permutations(a, t::Integer)
-    if t < 0
-        t = length(a) + 1
+    if t == 0
+        # Correct behavior for a permutation of length 0 is a vector containing a single empty vector 
+        return [Vector{eltype(a)}()]
+    elseif t == 1
+        # Easy case, just return each element in its own vector
+        return [[ai] for ai in a]
+    elseif t < 0 || t > length(a)
+        # Correct behavior for a permutation of these lengths is a an empty vector (of the correct type)
+        return Vector{Vector{eltype(a)}}()
     end
-    Permutations(a, t)
+    return Permutations(a, t)
 end
 
-function Base.iterate(p::Permutations, s = collect(1:length(p.a)))
-    (!isempty(s) && max(s[1], p.t) > length(p.a) || (isempty(s) && p.t > 0)) && return
-    nextpermutation(p.a, p.t ,s)
-end
 
 function nextpermutation(m, t, state)
     perm = [m[state[i]] for i in 1:t]
     n = length(state)
     if t <= 0
-        return(perm, [n+1])
+        return (perm, [n + 1])
     end
     s = copy(state)
     if t < n
         j = t + 1
-        while j <= n &&  s[t] >= s[j]; j+=1; end
+        while j <= n && s[t] >= s[j]
+            j += 1
+        end
     end
     if t < n && j <= n
         s[t], s[j] = s[j], s[t]
     else
         if t < n
-            reverse!(s, t+1)
+            reverse!(s, t + 1)
         end
         i = t - 1
-        while i>=1 && s[i] >= s[i+1]; i -= 1; end
+        while i >= 1 && s[i] >= s[i+1]
+            i -= 1
+        end
         if i > 0
             j = n
-            while j>i && s[i] >= s[j]; j -= 1; end
+            while j > i && s[i] >= s[j]
+                j -= 1
+            end
             s[i], s[j] = s[j], s[i]
-            reverse!(s, i+1)
+            reverse!(s, i + 1)
         else
-            s[1] = n+1
+            s[1] = n + 1
         end
     end
     return (perm, s)
@@ -94,18 +150,18 @@ function Base.length(c::MultiSetPermutations)
     else
         g = [factorial(i) for i in 0:t]
     end
-    p = [g[t+1]; zeros(Float64,t)]
+    p = [g[t+1]; zeros(Float64, t)]
     for i in 1:length(c.f)
         f = c.f[i]
         if i == 1
             for j in 1:min(f, t)
-                p[j+1] = g[t+1]/g[j+1]
+                p[j+1] = g[t+1] / g[j+1]
             end
         else
             for j in t:-1:1
                 q = 0
-                for k in (j+1):-1:max(1,j+1-f)
-                    q += p[k]/g[j+2-k]
+                for k in (j+1):-1:max(1, j + 1 - f)
+                    q += p[k] / g[j+2-k]
                 end
                 p[j+1] = q
             end
@@ -134,7 +190,7 @@ function multiset_permutations(m, f::Vector{<:Integer}, t::Integer)
     MultiSetPermutations(m, f, t, ref)
 end
 
-function Base.iterate(p::MultiSetPermutations, s = p.ref)
+function Base.iterate(p::MultiSetPermutations, s=p.ref)
     (!isempty(s) && max(s[1], p.t) > length(p.ref) || (isempty(s) && p.t > 0)) && return
     nextpermutation(p.m, p.t, s)
 end
@@ -151,7 +207,7 @@ function nthperm!(a::AbstractVector, k::Integer)
     f = factorial(oftype(k, n))
     0 < k <= f || throw(ArgumentError("permutation k must satisfy 0 < k ≤ $f, got $k"))
     k -= 1 # make k 1-indexed
-    for i=1:n-1
+    for i = 1:n-1
         f ÷= n - i + 1
         j = k ÷ f
         k -= j * f
@@ -182,7 +238,7 @@ function nthperm(p::AbstractVector{<:Integer})
     isperm(p) || throw(ArgumentError("argument is not a permutation"))
     k, n = 1, length(p)
     for i = 1:n-1
-        f = factorial(n-i)
+        f = factorial(n - i)
         for j = i+1:n
             k += ifelse(p[j] < p[i], f, 0)
         end
@@ -193,9 +249,10 @@ end
 
 # Parity of permutations
 
-const levicivita_lut = cat([0 0  0;  0 0 1; 0 -1 0],
-                           [0 0 -1;  0 0 0; 1  0 0],
-                           [0 1  0; -1 0 0; 0  0 0]; dims=3)
+const levicivita_lut = cat([0 0 0; 0 0 1; 0 -1 0],
+                           [0 0 -1; 0 0 0; 1 0 0],
+                           [0 1 0; -1 0 0; 0 0 0];
+                           dims=3)
 
 """
     levicivita(p)

test/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/combinations.jl

@@ -1,43 +1,46 @@
-@test [combinations([])...] == []
-@test [combinations(['a', 'b', 'c'])...] == [['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
-
-@test [combinations("abc",3)...] == [['a','b','c']]
-@test [combinations("abc",2)...] == [['a','b'],['a','c'],['b','c']]
-@test [combinations("abc",1)...] == [['a'],['b'],['c']]
-@test [combinations("abc",0)...] == [[]]
-@test [combinations("abc",-1)...] == []
-
-@test filter(x->iseven(x[1]),[combinations([1,2,3],2)...]) == Any[[2,3]]
-
-# multiset_combinations
-@test [multiset_combinations("aabc", 5)...] == Any[]
-@test [multiset_combinations("aabc", 2)...] == Any[['a','a'],['a','b'],['a','c'],['b','c']]
-@test [multiset_combinations("aabc", 1)...] == Any[['a'],['b'],['c']]
-@test [multiset_combinations("aabc", 0)...] == Any[Char[]]
-@test [multiset_combinations("aabc", -1)...] == Any[]
-@test [multiset_combinations("", 1)...] == Any[]
-@test [multiset_combinations("", 0)...] == Any[Char[]]
-@test [multiset_combinations("", -1)...] == Any[]
-
-# with_replacement_combinations
-@test [with_replacement_combinations("abc", 2)...] == Any[['a','a'],['a','b'],['a','c'],
-                                                              ['b','b'],['b','c'],['c','c']]
-@test [with_replacement_combinations("abc", 1)...] == Any[['a'],['b'],['c']]
-@test [with_replacement_combinations("abc", 0)...] == Any[Char[]]
-@test [with_replacement_combinations("abc", -1)...] == Any[]
-@test [with_replacement_combinations("", 1)...] == Any[]
-@test [with_replacement_combinations("", 0)...] == Any[Char[]]
-@test [with_replacement_combinations("", -1)...] == Any[]
-
-
-#cool-lex iterator
-@test_throws DomainError [CoolLexCombinations(-1, 1)...]
-@test_throws DomainError [CoolLexCombinations(5, 0)...]
-@test [CoolLexCombinations(4,2)...] == Vector[[1,2], [2,3], [1,3], [2,4], [3,4], [1,4]]
-@test isa(iterate(CoolLexCombinations(1000, 20))[2], Combinatorics.CoolLexIterState{BigInt})
-
-# Power set
-@test collect(powerset([])) == Any[[]]
-@test collect(powerset(['a', 'b', 'c'])) == Any[[],['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
-@test collect(powerset(['a', 'b', 'c'], 1)) == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
-@test collect(powerset(['a', 'b', 'c'], 1, 2)) == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c']]
+@testset "combinations" begin
+    @test [combinations([])...] == []
+    @test [combinations(['a', 'b', 'c'])...] == [['a'], ['b'], ['c'], ['a', 'b'], ['a', 'c'], ['b', 'c'], ['a', 'b', 'c']]
+
+    @test [combinations("abc", 3)...] == [['a', 'b', 'c']]
+    @test [combinations("abc", 2)...] == [['a', 'b'], ['a', 'c'], ['b', 'c']]
+    @test [combinations("abc", 1)...] == [['a'], ['b'], ['c']]
+    @test [combinations("abc", 0)...] == [[]]
+    @test [combinations("abc", -1)...] == []
+
+    @test filter(x -> iseven(x[1]), [combinations([1, 2, 3], 2)...]) == Any[[2, 3]]
+
+    # multiset_combinations
+    @test [multiset_combinations("aabc", 5)...] == Any[]
+    @test [multiset_combinations("aabc", 2)...] == Any[['a', 'a'], ['a', 'b'], ['a', 'c'], ['b', 'c']]
+    @test [multiset_combinations("aabc", 1)...] == Any[['a'], ['b'], ['c']]
+    @test [multiset_combinations("aabc", 0)...] == Any[Char[]]
+    @test [multiset_combinations("aabc", -1)...] == Any[]
+    @test [multiset_combinations("", 1)...] == Any[]
+    @test [multiset_combinations("", 0)...] == Any[Char[]]
+    @test [multiset_combinations("", -1)...] == Any[]
+
+    # with_replacement_combinations
+    @test [with_replacement_combinations("abc", 2)...] == Any[['a', 'a'], ['a', 'b'], ['a', 'c'],
+        ['b', 'b'], ['b', 'c'], ['c', 'c']]
+    @test [with_replacement_combinations("abc", 1)...] == Any[['a'], ['b'], ['c']]
+    @test [with_replacement_combinations("abc", 0)...] == Any[Char[]]
+    @test [with_replacement_combinations("abc", -1)...] == Any[]
+    @test [with_replacement_combinations("", 1)...] == Any[]
+    @test [with_replacement_combinations("", 0)...] == Any[Char[]]
+    @test [with_replacement_combinations("", -1)...] == Any[]
+
+
+    #cool-lex iterator
+    @test_throws DomainError [CoolLexCombinations(-1, 1)...]
+    @test_throws DomainError [CoolLexCombinations(5, 0)...]
+    @test [CoolLexCombinations(4, 2)...] == Vector[[1, 2], [2, 3], [1, 3], [2, 4], [3, 4], [1, 4]]
+    @test isa(iterate(CoolLexCombinations(1000, 20))[2], Combinatorics.CoolLexIterState{BigInt})
+
+    # Power set
+    @test collect(powerset([])) == Any[[]]
+    @test collect(powerset(['a', 'b', 'c'])) == Any[[], ['a'], ['b'], ['c'], ['a', 'b'], ['a', 'c'], ['b', 'c'], ['a', 'b', 'c']]
+    @test collect(powerset(['a', 'b', 'c'], 1)) == Any[['a'], ['b'], ['c'], ['a', 'b'], ['a', 'c'], ['b', 'c'], ['a', 'b', 'c']]
+    @test collect(powerset(['a', 'b', 'c'], 1, 2)) == Any[['a'], ['b'], ['c'], ['a', 'b'], ['a', 'c'], ['b', 'c']]
+
+end

test/factorials.jl

@@ -1,38 +1,41 @@
-@test factorial(7,3) == 7*6*5*4
-@test_throws DomainError factorial(3,7)
-@test_throws DomainError factorial(-3,-7)
-@test_throws DomainError factorial(-7,-3)
-#JuliaLang/julia#9943
-@test factorial(big(100), (80)) == 1303995018204712451095685346159820800000
-#JuliaLang/julia#9950
-@test_throws OverflowError factorial(1000,80)
-
-# derangement
-@test derangement(4) == subfactorial(4) == 9
-@test derangement(24) == parse(BigInt,"228250211305338670494289")
-
-# partialderangement
-@test partialderangement(7, 3) == 315
-@test_throws DomainError partialderangement(8, 9)
-@test_throws DomainError partialderangement(-8, 0)
-
-# doublefactorial
-@test doublefactorial(70) == parse(BigInt,"355044260642859198243475901411974413130137600000000")
-@test_throws DomainError doublefactorial(-1)
-
-# hyperfactorial
-@test hyperfactorial(8) == parse(BigInt,"55696437941726556979200000")
-@test hyperfactorial(0) == parse(BigInt,"1")
-@test hyperfactorial(1) == parse(BigInt,"1")
-@test hyperfactorial(2) == parse(BigInt,"4")
-
-# multifactorial
-@test multifactorial(40,2) == doublefactorial(40)
-@test_throws DomainError multifactorial(-1,1)
-
-# multinomial
-@test multinomial(1,4,4,2) == 34650
-
-# primorial
-@test primorial(17) == 510510
-@test_throws DomainError primorial(-1)
+@testset "factorials" begin
+    @test factorial(7, 3) == 7 * 6 * 5 * 4
+    @test_throws DomainError factorial(3, 7)
+    @test_throws DomainError factorial(-3, -7)
+    @test_throws DomainError factorial(-7, -3)
+    #JuliaLang/julia#9943
+    @test factorial(big(100), (80)) == 1303995018204712451095685346159820800000
+    #JuliaLang/julia#9950
+    @test_throws OverflowError factorial(1000, 80)
+
+    # derangement
+    @test derangement(4) == subfactorial(4) == 9
+    @test derangement(24) == parse(BigInt, "228250211305338670494289")
+
+    # partialderangement
+    @test partialderangement(7, 3) == 315
+    @test_throws DomainError partialderangement(8, 9)
+    @test_throws DomainError partialderangement(-8, 0)
+
+    # doublefactorial
+    @test doublefactorial(70) == parse(BigInt, "355044260642859198243475901411974413130137600000000")
+    @test_throws DomainError doublefactorial(-1)
+
+    # hyperfactorial
+    @test hyperfactorial(8) == parse(BigInt, "55696437941726556979200000")
+    @test hyperfactorial(0) == parse(BigInt, "1")
+    @test hyperfactorial(1) == parse(BigInt, "1")
+    @test hyperfactorial(2) == parse(BigInt, "4")
+
+    # multifactorial
+    @test multifactorial(40, 2) == doublefactorial(40)
+    @test_throws DomainError multifactorial(-1, 1)
+
+    # multinomial
+    @test multinomial(1, 4, 4, 2) == 34650
+
+    # primorial
+    @test primorial(17) == 510510
+    @test_throws DomainError primorial(-1)
+
+end