diff --git a/README.md b/README.md
index da278fc..2965b53 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,5 @@
# Combinatorics
-[](http://pkg.julialang.org/?pkg=Combinatorics&ver=0.3)
[](http://pkg.julialang.org/?pkg=Combinatorics&ver=0.4)
[](https://travis-ci.org/JuliaLang/Combinatorics.jl)
[](https://coveralls.io/github/JuliaLang/Combinatorics.jl?branch=master)
@@ -11,18 +10,18 @@ combinatorics and permutations. As overflows are expected even for low values,
most of the functions always return `BigInt`, and are marked as such below.
This library provides the following functions:
- - `bell(n)`: returns the n-th Bell number; always returns a `BigInt`;
- - `catalan(n)`: returns the n-th Catalan number; always returns a `BigInt`;
- - `combinations(a)`: returns combinations of all order by chaining calls to `Base.combinations(a,n);
+ - `bellnum(n)`: returns the n-th Bell number; always returns a `BigInt`;
+ - `catalannum(n)`: returns the n-th Catalan number; always returns a `BigInt`;
+ - `combinations(a)`: returns combinations of all order by chaining calls to `Base.combinations(a,n);
- `derangement(n)`/`subfactorial(n)`: returns the number of permutations of n with no fixed points; always returns a `BigInt`;
- `doublefactorial(n)`: returns the double factorial n!!; always returns a `BigInt`;
- `fibonacci(n)`: the n-th Fibonacci number; always returns a `BigInt`;
- `hyperfactorial(n)`: the n-th hyperfactorial, i.e. prod([i^i for i = 2:n]; always returns a `BigInt`;
- `integer_partitions(n)`: returns a `Vector{Int}` consisting of the partitions of the number `n`.
- `jacobisymbol(a,b)`: returns the Jacobi symbol (a/b);
- - `lassalle(n)`: returns the nth Lassalle number An defined in [arXiv:1009.4225](http://arxiv.org/abs/1009.4225) ([OEIS A180874](http://oeis.org/A180874)); always returns a `BigInt`;
+ - `lassallenum(n)`: returns the nth Lassalle number An defined in [arXiv:1009.4225](http://arxiv.org/abs/1009.4225) ([OEIS A180874](http://oeis.org/A180874)); always returns a `BigInt`;
- `legendresymbol(a,p)`: returns the Legendre symbol (a/p);
- - `lucas(n)`: the n-th Lucas number; always returns a `BigInt`;
+ - `lucasnum(n)`: the n-th Lucas number; always returns a `BigInt`;
- `multifactorial(n)`: returns the m-multifactorial n(!^m); always returns a `BigInt`;
- `multinomial(k...)`: receives a tuple of `k_1, ..., k_n` and calculates the multinomial coefficient `(n k)`, where `n = sum(k)`; returns a `BigInt` only if given a `BigInt`;
- `primorial(n)`: returns the product of all positive prime numbers <= n; always returns a `BigInt`;
diff --git a/REQUIRE b/REQUIRE
index fc421e0..77474e8 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,4 +1,3 @@
-julia 0.4
-Compat
+julia 0.5-
Polynomials
Iterators
diff --git a/src/Combinatorics.jl b/src/Combinatorics.jl
index b9a1df8..44fea2b 100644
--- a/src/Combinatorics.jl
+++ b/src/Combinatorics.jl
@@ -2,158 +2,13 @@ module Combinatorics
using Compat, Polynomials, Iterators
-import Base:combinations
-
-export bell,
- derangement,
- doublefactorial,
- fibonacci,
- hyperfactorial,
- jacobisymbol,
- lassalle,
- legendresymbol,
- lucas,
- multifactorial,
- multinomial,
- primorial,
- stirlings1,
- subfactorial
+import Base: start, next, done, length, eltype
+include("numbers.jl")
+include("factorials.jl")
+include("combinations.jl")
+include("permutations.jl")
include("partitions.jl")
include("youngdiagrams.jl")
-# Returns the n-th Bell number
-function bell(bn::Integer)
- if bn < 0
- throw(DomainError())
- else
- n = BigInt(bn)
- end
- list = Array(BigInt, div(n*(n+1), 2))
- list[1] = 1
- for i = 2:n
- beg = div(i*(i-1),2)
- list[beg+1] = list[beg]
- for j = 2:i
- list[beg+j] = list[beg+j-1]+list[beg+j-i]
- end
- end
- return list[end]
-end
-
-# Returns the n-th Catalan number
-function catalan(bn::Integer)
- if bn<0
- throw(DomainError())
- else
- n = BigInt(bn)
- end
- div(binomial(2*n, n), (n + 1))
-end
-
-#generate combinations of all orders, chaining of order iterators is eager,
-#but sequence at each order is lazy
-combinations(a) = chain([combinations(a,k) for k=1:length(a)]...)
-
-# The number of permutations of n with no fixed points (subfactorial)
-function derangement(sn::Integer)
- n = BigInt(sn)
- return num(factorial(n)*sum([(-1)^k//factorial(k) for k=0:n]))
-end
-subfactorial(n::Integer) = derangement(n)
-
-function doublefactorial(n::Integer)
- if n < 0
- throw(DomainError())
- end
- z = BigInt()
- ccall((:__gmpz_2fac_ui, :libgmp), Void,
- (Ptr{BigInt}, UInt), &z, @compat(UInt(n)))
- return z
-end
-
-function fibonacci(n::Integer)
- if n < 0
- throw(DomainError())
- end
- z = BigInt()
- ccall((:__gmpz_fib_ui, :libgmp), Void,
- (Ptr{BigInt}, UInt), &z, @compat(UInt(n)))
- return z
-end
-
-# Hyperfactorial
-hyperfactorial(n::Integer) = prod([i^i for i = BigInt(2):n])
-
-function jacobisymbol(a::Integer, b::Integer)
- ba = BigInt(a)
- bb = BigInt(b)
- return ccall((:__gmpz_jacobi, :libgmp), Cint,
- (Ptr{BigInt}, Ptr{BigInt}), &ba, &bb)
-end
-
-#Computes Lassalle's sequence
-#OEIS entry A180874
-function lassalle(m::Integer)
- A = ones(BigInt,m)
- for n=2:m
- A[n]=(-1)^(n-1) * (catalan(n) + sum([(-1)^j*binomial(2n-1, 2j-1)*A[j]*catalan(n-j) for j=1:n-1]))
- end
- A[m]
-end
-
-function legendresymbol(a::Integer, b::Integer)
- ba = BigInt(a)
- bb = BigInt(b)
- return ccall((:__gmpz_legendre, :libgmp), Cint,
- (Ptr{BigInt}, Ptr{BigInt}), &ba, &bb)
-end
-
-function lucas(n::Integer)
- if n < 0
- throw(DomainError())
- end
- z = BigInt()
- ccall((:__gmpz_lucnum_ui, :libgmp), Void,
- (Ptr{BigInt}, UInt), &z, @compat(UInt(n)))
- return z
-end
-
-function multifactorial(n::Integer, m::Integer)
- if n < 0
- throw(DomainError())
- end
- z = BigInt()
- ccall((:__gmpz_mfac_uiui, :libgmp), Void,
- (Ptr{BigInt}, UInt, UInt), &z, @compat(UInt(n)), @compat(UInt(m)))
- return z
-end
-
-# Multinomial coefficient where n = sum(k)
-function multinomial(k...)
- s = 0
- result = 1
- for i in k
- s += i
- result *= binomial(s, i)
- end
- result
-end
-
-function primorial(n::Integer)
- if n < 0
- throw(DomainError())
- end
- z = BigInt()
- ccall((:__gmpz_primorial_ui, :libgmp), Void,
- (Ptr{BigInt}, UInt), &z, @compat(UInt(n)))
- return z
-end
-
-# Returns s(n, k), the signed Stirling number of first kind
-function stirlings1(n::Integer, k::Integer)
- p = poly(0:(n-1))
- p[n - k + 1]
-end
-
-end # module
+end #module
diff --git a/src/combinations.jl b/src/combinations.jl
new file mode 100644
index 0000000..ca117cc
--- /dev/null
+++ b/src/combinations.jl
@@ -0,0 +1,52 @@
+export combinations
+
+#The Combinations iterator
+
+immutable Combinations{T}
+ a::T
+ t::Int
+end
+
+start(c::Combinations) = [1:c.t;]
+function next(c::Combinations, s)
+ comb = [c.a[si] for si in s]
+ if c.t == 0
+ # special case to generate 1 result for t==0
+ return (comb,[length(c.a)+2])
+ end
+ s = copy(s)
+ for i = length(s):-1:1
+ s[i] += 1
+ if s[i] > (length(c.a) - (length(s)-i))
+ continue
+ end
+ for j = i+1:endof(s)
+ s[j] = s[j-1]+1
+ end
+ break
+ end
+ (comb,s)
+end
+done(c::Combinations, s) = !isempty(s) && s[1] > length(c.a)-c.t+1
+
+length(c::Combinations) = binomial(length(c.a),c.t)
+
+eltype{T}(::Type{Combinations{T}}) = Vector{eltype(T)}
+
+"Generate all combinations of `n` elements from an indexable object. Because the number of combinations can be very large, this function returns an iterator object. Use `collect(combinations(array,n))` to get an array of all combinations.
+"
+function combinations(a, t::Integer)
+ if t < 0
+ # generate 0 combinations for negative argument
+ t = length(a)+1
+ end
+ Combinations(a, t)
+end
+
+
+"""
+generate combinations of all orders, chaining of order iterators is eager,
+but sequence at each order is lazy
+"""
+combinations(a) = chain([combinations(a,k) for k=1:length(a)]...)
+
diff --git a/src/factorials.jl b/src/factorials.jl
new file mode 100644
index 0000000..312835f
--- /dev/null
+++ b/src/factorials.jl
@@ -0,0 +1,82 @@
+#Factorials and elementary coefficients
+
+export
+ derangement,
+ factorial,
+ subfactorial,
+ doublefactorial,
+ hyperfactorial,
+ multifactorial,
+ gamma,
+ primorial,
+ multinomial
+
+import Base: factorial
+
+"computes n!/k!"
+function factorial{T<:Integer}(n::T, k::T)
+ if k < 0 || n < 0 || k > n
+ throw(DomainError())
+ end
+ f = one(T)
+ while n > k
+ f = Base.checked_mul(f,n)
+ n -= 1
+ end
+ return f
+end
+factorial(n::Integer, k::Integer) = factorial(promote(n, k)...)
+
+
+"The number of permutations of n with no fixed points (subfactorial)"
+function derangement(sn::Integer)
+ n = BigInt(sn)
+ return num(factorial(n)*sum([(-1)^k//factorial(k) for k=0:n]))
+end
+subfactorial(n::Integer) = derangement(n)
+
+function doublefactorial(n::Integer)
+ if n < 0
+ throw(DomainError())
+ end
+ z = BigInt()
+ ccall((:__gmpz_2fac_ui, :libgmp), Void,
+ (Ptr{BigInt}, UInt), &z, UInt(n))
+ return z
+end
+
+# Hyperfactorial
+hyperfactorial(n::Integer) = prod([i^i for i = BigInt(2):n])
+
+function multifactorial(n::Integer, m::Integer)
+ if n < 0
+ throw(DomainError())
+ end
+ z = BigInt()
+ ccall((:__gmpz_mfac_uiui, :libgmp), Void,
+ (Ptr{BigInt}, UInt, UInt), &z, UInt(n), UInt(m))
+ return z
+end
+
+function primorial(n::Integer)
+ if n < 0
+ throw(DomainError())
+ end
+ z = BigInt()
+ ccall((:__gmpz_primorial_ui, :libgmp), Void,
+ (Ptr{BigInt}, UInt), &z, UInt(n))
+ return z
+end
+
+"Multinomial coefficient where n = sum(k)"
+function multinomial(k...)
+ s = 0
+ result = 1
+ @inbounds for i in k
+ s += i
+ result *= binomial(s, i)
+ end
+ result
+end
+
+
diff --git a/src/numbers.jl b/src/numbers.jl
new file mode 100644
index 0000000..354bb83
--- /dev/null
+++ b/src/numbers.jl
@@ -0,0 +1,93 @@
+#Special named numbers and symbols
+
+export bellnum,
+ catalannum,
+ fibonaccinum,
+ jacobisymbol,
+ lassallenum,
+ legendresymbol,
+ lucasnum,
+ stirlings1
+
+"Returns the n-th Bell number"
+function bellnum(bn::Integer)
+ if bn < 0
+ throw(DomainError())
+ else
+ n = BigInt(bn)
+ end
+ list = Array(BigInt, div(n*(n+1), 2))
+ list[1] = 1
+ for i = 2:n
+ beg = div(i*(i-1),2)
+ list[beg+1] = list[beg]
+ for j = 2:i
+ list[beg+j] = list[beg+j-1]+list[beg+j-i]
+ end
+ end
+ return list[end]
+end
+
+"Returns the n-th Catalan number"
+function catalannum(bn::Integer)
+ if bn<0
+ throw(DomainError())
+ else
+ n = BigInt(bn)
+ end
+ div(binomial(2*n, n), (n + 1))
+end
+
+function fibonaccinum(n::Integer)
+ if n < 0
+ throw(DomainError())
+ end
+ z = BigInt()
+ ccall((:__gmpz_fib_ui, :libgmp), Void,
+ (Ptr{BigInt}, UInt), &z, UInt(n))
+ return z
+end
+
+
+function jacobisymbol(a::Integer, b::Integer)
+ ba = BigInt(a)
+ bb = BigInt(b)
+ return ccall((:__gmpz_jacobi, :libgmp), Cint,
+ (Ptr{BigInt}, Ptr{BigInt}), &ba, &bb)
+end
+
+"""
+Computes Lassalle's sequence
+OEIS entry A180874
+"""
+function lassallenum(m::Integer)
+ A = ones(BigInt,m)
+ for n=2:m
+ A[n]=(-1)^(n-1) * (catalannum(n) + sum([(-1)^j*binomial(2n-1, 2j-1)*A[j]*catalannum(n-j) for j=1:n-1]))
+ end
+ A[m]
+end
+
+function legendresymbol(a::Integer, b::Integer)
+ ba = BigInt(a)
+ bb = BigInt(b)
+ return ccall((:__gmpz_legendre, :libgmp), Cint,
+ (Ptr{BigInt}, Ptr{BigInt}), &ba, &bb)
+end
+
+function lucasnum(n::Integer)
+ if n < 0
+ throw(DomainError())
+ end
+ z = BigInt()
+ ccall((:__gmpz_lucnum_ui, :libgmp), Void,
+ (Ptr{BigInt}, UInt), &z, UInt(n))
+ return z
+end
+
+"Returns s(n, k), the signed Stirling number of first kind"
+function stirlings1(n::Integer, k::Integer)
+ p = poly(0:(n-1))
+ p[n - k + 1]
+end
+
diff --git a/src/partitions.jl b/src/partitions.jl
index 48ccd1b..8d5b60e 100644
--- a/src/partitions.jl
+++ b/src/partitions.jl
@@ -1,8 +1,396 @@
-#Generative algorithms
+#Partitions
-export cool_lex, integer_partitions, ncpartitions
+export
+ cool_lex,
+ integer_partitions,
+ ncpartitions,
+ partitions,
+ prevprod
+ #nextprod,
-# Lists the partitions of the number n, the order is consistent with GAP
+
+#integer partitions
+
+immutable IntegerPartitions
+ n::Int
+end
+
+length(p::IntegerPartitions) = npartitions(p.n)
+
+start(p::IntegerPartitions) = Int[]
+done(p::IntegerPartitions, xs) = length(xs) == p.n
+next(p::IntegerPartitions, xs) = (xs = nextpartition(p.n,xs); (xs,xs))
+
+"""
+Generate all integer arrays that sum to `n`. Because the number of partitions can be very large, this function returns an iterator object. Use `collect(partitions(n))` to get an array of all partitions. The number of partitions to generate can be efficiently computed using `length(partitions(n))`.
+"""
+partitions(n::Integer) = IntegerPartitions(n)
+
+
+
+function nextpartition(n, as)
+ if isempty(as); return Int[n]; end
+
+ xs = similar(as,0)
+ sizehint!(xs,length(as)+1)
+
+ for i = 1:length(as)-1
+ if as[i+1] == 1
+ x = as[i]-1
+ push!(xs, x)
+ n -= x
+ while n > x
+ push!(xs, x)
+ n -= x
+ end
+ push!(xs, n)
+
+ return xs
+ end
+ push!(xs, as[i])
+ n -= as[i]
+ end
+ push!(xs, as[end]-1)
+ push!(xs, 1)
+
+ xs
+end
+
+let _npartitions = Dict{Int,Int}()
+ global npartitions
+ function npartitions(n::Int)
+ if n < 0
+ 0
+ elseif n < 2
+ 1
+ elseif (np = get(_npartitions, n, 0)) > 0
+ np
+ else
+ np = 0
+ sgn = 1
+ for k = 1:n
+ np += sgn * (npartitions(n-k*(3k-1)>>1) + npartitions(n-k*(3k+1)>>1))
+ sgn = -sgn
+ end
+ _npartitions[n] = np
+ end
+ end
+end
+
+# Algorithm H from TAoCP 7.2.1.4
+# Partition n into m parts
+# in colex order (lexicographic by reflected sequence)
+
+immutable FixedPartitions
+ n::Int
+ m::Int
+end
+
+length(f::FixedPartitions) = npartitions(f.n,f.m)
+
+"""
+Generate all arrays of `m` integers that sum to `n`. Because the number of partitions can be very large, this function returns an iterator object. Use `collect(partitions(n,m))` to get an array of all partitions. The number of partitions to generate can be efficiently computed using `length(partitions(n,m))`.
+"""
+partitions(n::Integer, m::Integer) = n >= 1 && m >= 1 ? FixedPartitions(n,m) : throw(DomainError())
+
+start(f::FixedPartitions) = Int[]
+function done(f::FixedPartitions, s::Vector{Int})
+ f.m <= f.n || return true
+ isempty(s) && return false
+ return f.m == 1 || s[1]-1 <= s[end]
+end
+next(f::FixedPartitions, s::Vector{Int}) = (xs = nextfixedpartition(f.n,f.m,s); (xs,xs))
+
+function nextfixedpartition(n, m, bs)
+ as = copy(bs)
+ if isempty(as)
+ # First iteration
+ as = [n-m+1; ones(Int, m-1)]
+ elseif as[2] < as[1]-1
+ # Most common iteration
+ as[1] -= 1
+ as[2] += 1
+ else
+ # Iterate
+ local j
+ s = as[1]+as[2]-1
+ for j = 3:m
+ if as[j] < as[1]-1; break; end
+ s += as[j]
+ end
+ x = as[j] += 1
+ for k = j-1:-1:2
+ as[k] = x
+ s -= x
+ end
+ as[1] = s
+ end
+
+ return as
+end
+
+let _nipartitions = Dict{Tuple{Int,Int},Int}()
+ global npartitions
+ function npartitions(n::Int,m::Int)
+ if n < m || m == 0
+ 0
+ elseif n == m
+ 1
+ elseif (np = get(_nipartitions, (n,m), 0)) > 0
+ np
+ else
+ _nipartitions[(n,m)] = npartitions(n-1,m-1) + npartitions(n-m,m)
+ end
+ end
+end
+
+# Algorithm H from TAoCP 7.2.1.5
+# Set partitions
+
+immutable SetPartitions{T<:AbstractVector}
+ s::T
+end
+
+length(p::SetPartitions) = nsetpartitions(length(p.s))
+
+"""
+Generate all set partitions of the elements of an array, represented as arrays of arrays. Because the number of partitions can be very large, this function returns an iterator object. Use `collect(partitions(array))` to get an array of all partitions. The number of partitions to generate can be efficiently computed using `length(partitions(array))`.
+"""
+partitions(s::AbstractVector) = SetPartitions(s)
+
+start(p::SetPartitions) = (n = length(p.s); (zeros(Int32, n), ones(Int32, n-1), n, 1))
+done(p::SetPartitions, s) = s[1][1] > 0
+next(p::SetPartitions, s) = nextsetpartition(p.s, s...)
+
+function nextsetpartition(s::AbstractVector, a, b, n, m)
+ function makeparts(s, a, m)
+ temp = [ similar(s,0) for k = 0:m ]
+ for i = 1:n
+ push!(temp[a[i]+1], s[i])
+ end
+ filter!(x->!isempty(x), temp)
+ end
+
+ if isempty(s); return ([s], ([1], Int[], n, 1)); end
+
+ part = makeparts(s,a,m)
+
+ if a[end] != m
+ a[end] += 1
+ else
+ local j
+ for j = n-1:-1:1
+ if a[j] != b[j]
+ break
+ end
+ end
+ a[j] += 1
+ m = b[j] + (a[j] == b[j])
+ for k = j+1:n-1
+ a[k] = 0
+ b[k] = m
+ end
+ a[end] = 0
+ end
+
+ return (part, (a,b,n,m))
+
+end
+
+let _nsetpartitions = Dict{Int,Int}()
+ global nsetpartitions
+ function nsetpartitions(n::Int)
+ if n < 0
+ 0
+ elseif n < 2
+ 1
+ elseif (wn = get(_nsetpartitions, n, 0)) > 0
+ wn
+ else
+ wn = 0
+ for k = 0:n-1
+ wn += binomial(n-1,k)*nsetpartitions(n-1-k)
+ end
+ _nsetpartitions[n] = wn
+ end
+ end
+end
+
+immutable FixedSetPartitions{T<:AbstractVector}
+ s::T
+ m::Int
+end
+
+length(p::FixedSetPartitions) = nfixedsetpartitions(length(p.s),p.m)
+
+"""
+Generate all set partitions of the elements of an array into exactly m subsets, represented as arrays of arrays. Because the number of partitions can be very large, this function returns an iterator object. Use `collect(partitions(array,m))` to get an array of all partitions. The number of partitions into m subsets is equal to the Stirling number of the second kind and can be efficiently computed using `length(partitions(array,m))`.
+"""
+partitions(s::AbstractVector,m::Int) = length(s) >= 1 && m >= 1 ? FixedSetPartitions(s,m) : throw(DomainError())
+
+function start(p::FixedSetPartitions)
+ n = length(p.s)
+ m = p.m
+ m <= n ? (vcat(ones(Int, n-m),1:m), vcat(1,n-m+2:n), n) : (Int[], Int[], n)
+end
+# state consists of:
+# vector a of length n describing to which partition every element of s belongs
+# vector b of length n describing the first index b[i] that belongs to partition i
+# integer n
+
+done(p::FixedSetPartitions, s) = isempty(s[1]) || s[1][1] > 1
+next(p::FixedSetPartitions, s) = nextfixedsetpartition(p.s,p.m, s...)
+
+function nextfixedsetpartition(s::AbstractVector, m, a, b, n)
+ function makeparts(s, a)
+ part = [ similar(s,0) for k = 1:m ]
+ for i = 1:n
+ push!(part[a[i]], s[i])
+ end
+ return part
+ end
+
+ part = makeparts(s,a)
+
+ if m == 1
+ a[1] = 2
+ return (part, (a, b, n))
+ end
+
+ if a[end] != m
+ a[end] += 1
+ else
+ local j, k
+ for j = n-1:-1:1
+ if a[j]1
+ a[j]+=1
+ for p=j+1:n
+ if b[a[p]]!=p
+ a[p]=1
+ end
+ end
+ else
+ for k=m:-1:2
+ if b[k-1]= x
+#for integer n1, n2, n3
+#"""
+#function nextprod(a::Vector{Int}, x)
+# if x > typemax(Int)
+# throw(ArgumentError("unsafe for x > typemax(Int), got $x"))
+# end
+# k = length(a)
+# v = ones(Int, k) # current value of each counter
+# mx = [nextpow(ai,x) for ai in a] # maximum value of each counter
+# v[1] = mx[1] # start at first case that is >= x
+# p::widen(Int) = mx[1] # initial value of product in this case
+# best = p
+# icarry = 1
+#
+# while v[end] < mx[end]
+# if p >= x
+# best = p < best ? p : best # keep the best found yet
+# carrytest = true
+# while carrytest
+# p = div(p, v[icarry])
+# v[icarry] = 1
+# icarry += 1
+# p *= a[icarry]
+# v[icarry] *= a[icarry]
+# carrytest = v[icarry] > mx[icarry] && icarry < k
+# end
+# if p < x
+# icarry = 1
+# end
+# else
+# while p < x
+# p *= a[1]
+# v[1] *= a[1]
+# end
+# end
+# end
+# best = mx[end] < best ? mx[end] : best
+# return Int(best) # could overflow, but best to have predictable return type
+#end
+
+"""
+Previous integer not greater than `n` that can be written as $\prod k_i^{p_i}$ for integers $p_1$, $p_2$, etc.
+
+For a list of integers i1, i2, i3, find the largest
+ i1^n1 * i2^n2 * i3^n3 <= x
+for integer n1, n2, n3
+"""
+function prevprod(a::Vector{Int}, x)
+ if x > typemax(Int)
+ throw(ArgumentError("unsafe for x > typemax(Int), got $x"))
+ end
+ k = length(a)
+ v = ones(Int, k) # current value of each counter
+ mx = [nextpow(ai,x) for ai in a] # allow each counter to exceed p (sentinel)
+ first = Int(prevpow(a[1], x)) # start at best case in first factor
+ v[1] = first
+ p::widen(Int) = first
+ best = p
+ icarry = 1
+
+ while v[end] < mx[end]
+ while p <= x
+ best = p > best ? p : best
+ p *= a[1]
+ v[1] *= a[1]
+ end
+ if p > x
+ carrytest = true
+ while carrytest
+ p = div(p, v[icarry])
+ v[icarry] = 1
+ icarry += 1
+ p *= a[icarry]
+ v[icarry] *= a[icarry]
+ carrytest = v[icarry] > mx[icarry] && icarry < k
+ end
+ if p <= x
+ icarry = 1
+ end
+ end
+ end
+ best = x >= p > best ? p : best
+ return Int(best)
+end
+
+
+"Lists the partitions of the number n, the order is consistent with GAP"
function integer_partitions(n::Integer)
if n < 0
throw(DomainError())
@@ -24,20 +412,22 @@ function integer_partitions(n::Integer)
list
end
-# Produces (n,k)-combinations in cool-lex order
-#
-#Implements the cool-lex algorithm to generate (n,k)-combinations
-#@article{Ruskey:2009fk,
-# Author = {Frank Ruskey and Aaron Williams},
-# Doi = {10.1016/j.disc.2007.11.048},
-# Journal = {Discrete Mathematics},
-# Month = {September},
-# Number = {17},
-# Pages = {5305-5320},
-# Title = {The coolest way to generate combinations},
-# Url = {http://www.sciencedirect.com/science/article/pii/S0012365X07009570},
-# Volume = {309},
-# Year = {2009}}
+"""
+Produces (n,k)-combinations in cool-lex order
+
+Implements the cool-lex algorithm to generate (n,k)-combinations
+@article{Ruskey:2009fk,
+ Author = {Frank Ruskey and Aaron Williams},
+ Doi = {10.1016/j.disc.2007.11.048},
+ Journal = {Discrete Mathematics},
+ Month = {September},
+ Number = {17},
+ Pages = {5305-5320},
+ Title = {The coolest way to generate combinations},
+ Url = {http://www.sciencedirect.com/science/article/pii/S0012365X07009570},
+ Volume = {309},
+ Year = {2009}}
+"""
function cool_lex(n::Integer, t::Integer)
s = n-t
if n > 64 error("Not implemented for n > 64") end
diff --git a/src/permutations.jl b/src/permutations.jl
new file mode 100644
index 0000000..4a54ed5
--- /dev/null
+++ b/src/permutations.jl
@@ -0,0 +1,139 @@
+#Permutations
+
+export
+ levicivita,
+ nthperm!,
+ nthperm,
+ parity,
+ permutations
+
+#The basic permutations iterator
+
+immutable Permutations{T}
+ a::T
+end
+
+"""
+Generate all permutations of an indexable object. Because the number of permutations can be very large, this function returns an iterator object. Use `collect(permutations(array))` to get an array of all permutations.
+"""
+permutations(a) = Permutations(a)
+
+eltype{T}(::Type{Permutations{T}}) = Vector{eltype(T)}
+
+length(p::Permutations) = factorial(length(p.a))
+
+start(p::Permutations) = [1:length(p.a);]
+function next(p::Permutations, s)
+ perm = [p.a[si] for si in s]
+ if isempty(p.a)
+ # special case to generate 1 result for len==0
+ return (perm,[1])
+ end
+ s = copy(s)
+ k = length(s)-1
+ while k > 0 && s[k] > s[k+1]; k -= 1; end
+ if k == 0
+ s[1] = length(s)+1 # done
+ else
+ l = length(s)
+ while s[k] >= s[l]; l -= 1; end
+ s[k],s[l] = s[l],s[k]
+ reverse!(s,k+1)
+ end
+ (perm,s)
+end
+done(p::Permutations, s) = !isempty(s) && s[1] > length(p.a)
+
+
+"In-place version of nthperm."
+function nthperm!(a::AbstractVector, k::Integer)
+ k -= 1 # make k 1-indexed
+ k < 0 && throw(ArgumentError("permutation k must be ≥ 0, got $k"))
+ n = length(a)
+ n == 0 && return a
+ f = factorial(oftype(k, n-1))
+ for i=1:n-1
+ j = div(k, f) + 1
+ k = k % f
+ f = div(f, n-i)
+
+ j = j+i-1
+ elt = a[j]
+ for d = j:-1:i+1
+ a[d] = a[d-1]
+ end
+ a[i] = elt
+ end
+ a
+end
+
+"Compute the kth lexicographic permutation of the vector a."
+nthperm(a::AbstractVector, k::Integer) = nthperm!(copy(a),k)
+
+"Return the `k` that generated permutation `p`. Note that `nthperm(nthperm([1:n], k)) == k` for `1 <= k <= factorial(n)`."
+function nthperm{T<:Integer}(p::AbstractVector{T})
+ isperm(p) || throw(ArgumentError("argument is not a permutation"))
+ k, n = 1, length(p)
+ for i = 1:n-1
+ f = factorial(n-i)
+ for j = i+1:n
+ k += ifelse(p[j] < p[i], f, 0)
+ end
+ end
+ return k
+end
+
+
+# Parity of permutations
+
+const levicivita_lut = cat(3, [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])
+
+"""
+Levi-Civita symbol of a permutation.
+
+Returns 1 is the permutation is even, -1 if it is odd and 0 otherwise.
+
+The parity is computed by using the fact that a permutation is odd if and
+only if the number of even-length cycles is odd.
+"""
+function levicivita{T<:Integer}(p::AbstractVector{T})
+ n = length(p)
+
+ if n == 3
+ @inbounds valid = (0 < p[1] <= 3) * (0 < p[2] <= 3) * (0 < p[3] <= 3)
+ return valid ? levicivita_lut[p[1], p[2], p[3]] : 0
+ end
+
+ todo = trues(n)
+ first = 1
+ cycles = flips = 0
+
+ while cycles + flips < n
+ first = findnext(todo, first)
+ (todo[first] $= true) && return 0
+ j = p[first]
+ (0 < j <= n) || return 0
+ cycles += 1
+ while j ≠ first
+ (todo[j] $= true) && return 0
+ j = p[j]
+ (0 < j <= n) || return 0
+ flips += 1
+ end
+ end
+
+ return iseven(flips) ? 1 : -1
+end
+
+"""
+Computes the parity of a permutation using the levicivita function,
+so you can ask iseven(parity(p)). If p is not a permutation throws an error.
+"""
+function parity{T<:Integer}(p::AbstractVector{T})
+ epsilon = levicivita(p)
+ epsilon == 0 && throw(ArgumentError("Not a permutation"))
+ epsilon == 1 ? 0 : 1
+end
+
diff --git a/src/youngdiagrams.jl b/src/youngdiagrams.jl
index 1cc397d..e788ff6 100644
--- a/src/youngdiagrams.jl
+++ b/src/youngdiagrams.jl
@@ -98,7 +98,7 @@ end
# partition sequences #
#######################
-#Computes essential part of the partition sequence of lambda
+"Computes essential part of the partition sequence of lambda"
function partitionsequence(lambda::Partition)
Λ▔ = Int64[]
λ = [lambda; 0]
@@ -145,16 +145,18 @@ function MN1inner(R::Vector{Int64}, T::Dict, μ::Partition, t::Integer)
χ
end
-#Computes character $χ^λ(μ)$ of the partition μ in the λth irrep of the
-#symmetric group $S_n$
-#
-#Implements the Murnaghan-Nakayama algorithm as described in:
-# Dan Bernstein,
-# "The computational complexity of rules for the character table of Sn",
-# Journal of Symbolic Computation, vol. 37 iss. 6 (2004), pp 727-748.
-# doi:10.1016/j.jsc.2003.11.001
+"""
+Computes character $χ^λ(μ)$ of the partition μ in the λth irrep of the
+symmetric group $S_n$
+
+Implements the Murnaghan-Nakayama algorithm as described in:
+ Dan Bernstein,
+ "The computational complexity of rules for the character table of Sn",
+ Journal of Symbolic Computation, vol. 37 iss. 6 (2004), pp 727-748.
+ doi:10.1016/j.jsc.2003.11.001
+"""
function character(λ::Partition, μ::Partition)
- T = @compat Dict{Any,Any}(()=>0) #Sparse array implemented as dict
+ T = Dict(()=>0) #Sparse array implemented as dict
Λ▔ = partitionsequence(λ)
MN1inner(Λ▔, T, μ, 1)
end
diff --git a/test/basic.jl b/test/basic.jl
deleted file mode 100644
index c742390..0000000
--- a/test/basic.jl
+++ /dev/null
@@ -1,55 +0,0 @@
-using Combinatorics
-using Base.Test
-
-# catalan
-@test Combinatorics.catalan(5) == 42
-@test Combinatorics.catalan(30) == parse(BigInt,"3814986502092304")
-
-#combinations
-@test collect(combinations([])) == []
-@test collect(combinations(['a', 'b', 'c'])) == Vector{Char}[['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
-
-# derangement
-@test derangement(4) == 9
-@test derangement(24) == parse(BigInt,"228250211305338670494289")
-
-# doublefactorial
-@test doublefactorial(70) == parse(BigInt,"355044260642859198243475901411974413130137600000000")
-
-# fibonacci
-@test fibonacci(5) == 5
-@test fibonacci(101) == parse(BigInt,"573147844013817084101")
-
-# hyperfactorial
-@test hyperfactorial(8) == parse(BigInt,"55696437941726556979200000")
-
-# lassalle
-@test lassalle(14) == parse(BigInt,"270316008395632253340")
-
-# legendresymbol
-@test legendresymbol(1001,9907) == jacobisymbol(1001,9907) == -1
-
-# lucas
-@test lucas(10) == 123
-@test lucas(100) == parse(BigInt,"792070839848372253127")
-
-# multifactorial
-@test multifactorial(40,2) == doublefactorial(40)
-
-# multinomial
-@test multinomial(1,4,4,2) == 34650
-
-# primorial
-@test primorial(17) == 510510
-
-# stirlings1
-@test sum([abs(stirlings1(8, i)) for i = 0:8]) == factorial(8)
-
-# bell
-@test bell(7) == 877
-@test bell(42) == parse(BigInt,"35742549198872617291353508656626642567")
-@test_throws DomainError bell(-1)
-
-# integer_partitions
-@test integer_partitions(5) == Any[[1, 1, 1, 1, 1], [2, 1, 1, 1], [2, 2, 1], [3, 1, 1], [3, 2], [4, 1], [5]]
-@test_throws DomainError integer_partitions(-1)
diff --git a/test/combinations.jl b/test/combinations.jl
new file mode 100644
index 0000000..8ba59d7
--- /dev/null
+++ b/test/combinations.jl
@@ -0,0 +1,17 @@
+using Combinatorics
+using Base.Test
+
+import Combinatorics: combinations
+
+
+@test collect(combinations([])) == []
+@test collect(combinations(['a', 'b', 'c'])) == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
+
+@test collect(combinations("abc",3)) == Any[['a','b','c']]
+@test collect(combinations("abc",2)) == Any[['a','b'],['a','c'],['b','c']]
+@test collect(combinations("abc",1)) == Any[['a'],['b'],['c']]
+@test collect(combinations("abc",0)) == Any[Char[]]
+@test collect(combinations("abc",-1)) == Any[]
+
+@test collect(filter(x->(iseven(x[1])),combinations([1,2,3],2))) == Any[[2,3]]
+
diff --git a/test/factorials.jl b/test/factorials.jl
new file mode 100644
index 0000000..a9342bb
--- /dev/null
+++ b/test/factorials.jl
@@ -0,0 +1,32 @@
+using Combinatorics
+using Base.Test
+
+@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) == 9
+@test derangement(24) == parse(BigInt,"228250211305338670494289")
+
+# doublefactorial
+@test doublefactorial(70) == parse(BigInt,"355044260642859198243475901411974413130137600000000")
+
+# hyperfactorial
+@test hyperfactorial(8) == parse(BigInt,"55696437941726556979200000")
+
+# multifactorial
+@test multifactorial(40,2) == doublefactorial(40)
+
+# multinomial
+@test multinomial(1,4,4,2) == 34650
+
+# primorial
+@test primorial(17) == 510510
+
+
diff --git a/test/numbers.jl b/test/numbers.jl
new file mode 100644
index 0000000..84e0513
--- /dev/null
+++ b/test/numbers.jl
@@ -0,0 +1,30 @@
+using Combinatorics
+using Base.Test
+
+# catalan
+@test catalannum(5) == 42
+@test catalannum(30) == parse(BigInt,"3814986502092304")
+
+# fibonacci
+@test fibonaccinum(5) == 5
+@test fibonaccinum(101) == parse(BigInt,"573147844013817084101")
+
+# lassalle
+@test lassallenum(14) == parse(BigInt,"270316008395632253340")
+
+# legendresymbol
+@test legendresymbol(1001,9907) == jacobisymbol(1001,9907) == -1
+
+# lucas
+@test lucasnum(10) == 123
+@test lucasnum(100) == parse(BigInt,"792070839848372253127")
+
+# stirlings1
+@test sum([abs(stirlings1(8, i)) for i = 0:8]) == factorial(8)
+
+# bell
+@test bellnum(7) == 877
+@test bellnum(42) == parse(BigInt,"35742549198872617291353508656626642567")
+@test_throws DomainError bellnum(-1)
+
+
diff --git a/test/partitions.jl b/test/partitions.jl
new file mode 100644
index 0000000..7f90ddf
--- /dev/null
+++ b/test/partitions.jl
@@ -0,0 +1,28 @@
+using Combinatorics
+using Base.Test
+
+@test collect(partitions(4)) == Any[[4], [3,1], [2,2], [2,1,1], [1,1,1,1]]
+@test collect(partitions(8,3)) == Any[[6,1,1], [5,2,1], [4,3,1], [4,2,2], [3,3,2]]
+@test collect(partitions(8, 1)) == Any[[8]]
+@test collect(partitions(8, 9)) == []
+@test collect(partitions([1,2,3])) == Any[Any[[1,2,3]], Any[[1,2],[3]], Any[[1,3],[2]], Any[[1],[2,3]], Any[[1],[2],[3]]]
+@test collect(partitions([1,2,3,4],3)) == Any[Any[[1,2],[3],[4]], Any[[1,3],[2],[4]], Any[[1],[2,3],[4]],
+ Any[[1,4],[2],[3]], Any[[1],[2,4],[3]], Any[[1],[2],[3,4]]]
+@test collect(partitions([1,2,3,4],1)) == Any[Any[[1, 2, 3, 4]]]
+@test collect(partitions([1,2,3,4],5)) == []
+
+@test length(partitions(0)) == 1
+@test length(partitions(-1)) == 0
+@test length(collect(partitions(30))) == length(partitions(30))
+@test length(collect(partitions(90,4))) == length(partitions(90,4))
+@test length(collect(partitions('a':'h'))) == length(partitions('a':'h'))
+@test length(collect(partitions('a':'h',5))) == length(partitions('a':'h',5))
+
+# integer_partitions
+@test integer_partitions(5) == Any[[1, 1, 1, 1, 1], [2, 1, 1, 1], [2, 2, 1], [3, 1, 1], [3, 2], [4, 1], [5]]
+@test_throws DomainError integer_partitions(-1)
+
+@test_throws ArgumentError prevprod([2,3,5],Int128(typemax(Int))+1)
+@test prevprod([2,3,5],30) == 30
+@test prevprod([2,3,5],33) == 32
+
diff --git a/test/permutations.jl b/test/permutations.jl
new file mode 100644
index 0000000..e4275b1
--- /dev/null
+++ b/test/permutations.jl
@@ -0,0 +1,35 @@
+using Combinatorics
+using Base.Test
+
+import Combinatorics: levicivita, nthperm, nthperm!, parity, permutations
+
+@test collect(permutations("abc")) == Any[['a','b','c'],['a','c','b'],['b','a','c'],
+ ['b','c','a'],['c','a','b'],['c','b','a']]
+
+@test collect(filter(x->(iseven(x[1])),permutations([1,2,3]))) == Any[[2,1,3],[2,3,1]]
+@test collect(filter(x->(iseven(x[3])),permutations([1,2,3]))) == Any[[1,3,2],[3,1,2]]
+
+@test length(permutations(0)) == 1
+
+#nthperm!
+for n = 0:7, k = 1:factorial(n)
+ p = nthperm!([1:n;], k)
+ @test isperm(p)
+ @test nthperm(p) == k
+end
+#Euler Problem #24
+@test nthperm!([0:9;],1000000) == [2,7,8,3,9,1,5,4,6,0]
+
+@test nthperm([0,1,2],3) == [1,0,2]
+
+@test_throws ArgumentError parity([0])
+@test_throws ArgumentError parity([1,2,3,3])
+@test levicivita([1,1,2,3]) == 0
+@test levicivita([1]) == 1 && parity([1]) == 0
+@test map(levicivita, collect(permutations([1,2,3]))) == [1, -1, -1, 1, 1, -1]
+@test let p = [3, 4, 6, 10, 5, 2, 1, 7, 8, 9]; levicivita(p) == 1 && parity(p) == 0; end
+@test let p = [4, 3, 6, 10, 5, 2, 1, 7, 8, 9]; levicivita(p) == -1 && parity(p) == 1; end
+
+@test Combinatorics.nsetpartitions(-1) == 0
+@test collect(permutations([])) == Vector[[]]
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 9cb5de0..7acd6ef 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,10 @@
using Combinatorics
using Base.Test
-include("basic.jl")
+include("numbers.jl")
+include("factorials.jl")
+include("combinations.jl")
+include("permutations.jl")
+include("partitions.jl")
include("youngdiagrams.jl")