From 2e5614978c7b751c66a4a2d947344eb1ebf11db6 Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Fri, 6 Nov 2015 14:49:05 -0500 Subject: [PATCH 1/5] Major code reorganization - DROP 0.4 SUPPORT - Import most of base/combinatorics.jl (Ref: JuliaLang/julia#13897) - Move most of the special numbers to numbers.jl - Put combinations, permutations and partitions in their own files - Rename special numbers with ~num suffix. This renaming is particularly important for catalannum to avoid clashing with the Base.catalan irrational constant. --- README.md | 1 - REQUIRE | 2 +- src/Combinatorics.jl | 157 +------------------ src/combinations.jl | 51 ++++++ src/factorials.jl | 64 ++++++++ src/numbers.jl | 91 +++++++++++ src/partitions.jl | 362 ++++++++++++++++++++++++++++++++++++++++++- src/permutations.jl | 181 ++++++++++++++++++++++ test/basic.jl | 24 +-- 9 files changed, 766 insertions(+), 167 deletions(-) create mode 100644 src/combinations.jl create mode 100644 src/factorials.jl create mode 100644 src/numbers.jl create mode 100644 src/permutations.jl diff --git a/README.md b/README.md index da278fc..40ee390 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,5 @@ # Combinatorics -[![Combinatorics](http://pkg.julialang.org/badges/Combinatorics_0.3.svg)](http://pkg.julialang.org/?pkg=Combinatorics&ver=0.3) [![Combinatorics](http://pkg.julialang.org/badges/Combinatorics_0.4.svg)](http://pkg.julialang.org/?pkg=Combinatorics&ver=0.4) [![Build Status](https://travis-ci.org/JuliaLang/Combinatorics.jl.svg?branch=master)](https://travis-ci.org/JuliaLang/Combinatorics.jl) [![Coverage Status](https://coveralls.io/repos/JuliaLang/Combinatorics.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/JuliaLang/Combinatorics.jl?branch=master) diff --git a/REQUIRE b/REQUIRE index fc421e0..045ea20 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,4 +1,4 @@ -julia 0.4 +julia 0.5- Compat Polynomials Iterators diff --git a/src/Combinatorics.jl b/src/Combinatorics.jl index b9a1df8..be625fa 100644 --- a/src/Combinatorics.jl +++ b/src/Combinatorics.jl @@ -2,158 +2,11 @@ module Combinatorics using Compat, Polynomials, Iterators -import Base:combinations - -export bell, - derangement, - doublefactorial, - fibonacci, - hyperfactorial, - jacobisymbol, - lassalle, - legendresymbol, - lucas, - multifactorial, - multinomial, - primorial, - stirlings1, - subfactorial - +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..7a06e44 --- /dev/null +++ b/src/combinations.jl @@ -0,0 +1,51 @@ +export combinations + + +#The Combinations iterator +import Base: start, next, done, length + +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)} + + + +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..a4f2307 --- /dev/null +++ b/src/factorials.jl @@ -0,0 +1,64 @@ +#Factorials and elementary coefficients + +export + derangement, + subfactorial, + doublefactorial, + hyperfactorial, + multifactorial, + gamma, + primorial, + multinomial + +# 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 + +# 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, @compat(UInt(n)), @compat(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, @compat(UInt(n))) + 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 + + diff --git a/src/numbers.jl b/src/numbers.jl new file mode 100644 index 0000000..404e292 --- /dev/null +++ b/src/numbers.jl @@ -0,0 +1,91 @@ +#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, @compat(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, @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 + diff --git a/src/partitions.jl b/src/partitions.jl index 48ccd1b..641ca41 100644 --- a/src/partitions.jl +++ b/src/partitions.jl @@ -1,7 +1,367 @@ -#Generative algorithms +#Partitions export cool_lex, integer_partitions, ncpartitions + +#integer partitions + +immutable IntegerPartitions + n::Int +end + +length(p::IntegerPartitions) = npartitions(p.n) + +partitions(n::Integer) = IntegerPartitions(n) + +start(p::IntegerPartitions) = Int[] +done(p::IntegerPartitions, xs) = length(xs) == p.n +next(p::IntegerPartitions, xs) = (xs = nextpartition(p.n,xs); (xs,xs)) + +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) + +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)) + +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) + +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 + +# 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 diff --git a/src/permutations.jl b/src/permutations.jl new file mode 100644 index 0000000..badabd7 --- /dev/null +++ b/src/permutations.jl @@ -0,0 +1,181 @@ +immutable Permutations{T} + a::T +end + +eltype{T}(::Type{Permutations{T}}) = Vector{eltype(T)} + +length(p::Permutations) = factorial(length(p.a)) + +permutations(a) = Permutations(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) + +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 +nthperm(a::AbstractVector, k::Integer) = nthperm!(copy(a),k) + +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 + +function invperm(a::AbstractVector) + b = zero(a) # similar vector of zeros + n = length(a) + for i = 1:n + j = a[i] + ((1 <= j <= n) && b[j] == 0) || + throw(ArgumentError("argument is not a permutation")) + b[j] = i + end + b +end + +function isperm(A) + n = length(A) + used = falses(n) + for a in A + (0 < a <= n) && (used[a] $= true) || return false + end + true +end + +function permute!!{T<:Integer}(a, p::AbstractVector{T}) + count = 0 + start = 0 + while count < length(a) + ptr = start = findnext(p, start+1) + temp = a[start] + next = p[start] + count += 1 + while next != start + a[ptr] = a[next] + p[ptr] = 0 + ptr = next + next = p[next] + count += 1 + end + a[ptr] = temp + p[ptr] = 0 + end + a +end + +permute!(a, p::AbstractVector) = permute!!(a, copy!(similar(p), p)) + +function ipermute!!{T<:Integer}(a, p::AbstractVector{T}) + count = 0 + start = 0 + while count < length(a) + start = findnext(p, start+1) + temp = a[start] + next = p[start] + count += 1 + while next != start + temp_next = a[next] + a[next] = temp + temp = temp_next + ptr = p[next] + p[next] = 0 + next = ptr + count += 1 + end + a[next] = temp + p[next] = 0 + end + a +end + +ipermute!(a, p::AbstractVector) = ipermute!!(a, copy!(similar(p), p)) + + +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. +# 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. +# Returns 1 is the permutarion is even, -1 if it is odd and 0 otherwise. +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/test/basic.jl b/test/basic.jl index c742390..8be45e1 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -2,12 +2,12 @@ using Combinatorics using Base.Test # catalan -@test Combinatorics.catalan(5) == 42 -@test Combinatorics.catalan(30) == parse(BigInt,"3814986502092304") +@test catalannum(5) == 42 +@test catalannum(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']] +@test collect(Combinatorics.combinations([])) == [] +@test collect(Combinatorics.combinations(['a', 'b', 'c'])) == Vector{Char}[['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']] # derangement @test derangement(4) == 9 @@ -17,21 +17,21 @@ using Base.Test @test doublefactorial(70) == parse(BigInt,"355044260642859198243475901411974413130137600000000") # fibonacci -@test fibonacci(5) == 5 -@test fibonacci(101) == parse(BigInt,"573147844013817084101") +@test fibonaccinum(5) == 5 +@test fibonaccinum(101) == parse(BigInt,"573147844013817084101") # hyperfactorial @test hyperfactorial(8) == parse(BigInt,"55696437941726556979200000") # lassalle -@test lassalle(14) == parse(BigInt,"270316008395632253340") +@test lassallenum(14) == parse(BigInt,"270316008395632253340") # legendresymbol @test legendresymbol(1001,9907) == jacobisymbol(1001,9907) == -1 # lucas -@test lucas(10) == 123 -@test lucas(100) == parse(BigInt,"792070839848372253127") +@test lucasnum(10) == 123 +@test lucasnum(100) == parse(BigInt,"792070839848372253127") # multifactorial @test multifactorial(40,2) == doublefactorial(40) @@ -46,9 +46,9 @@ using Base.Test @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) +@test bellnum(7) == 877 +@test bellnum(42) == parse(BigInt,"35742549198872617291353508656626642567") +@test_throws DomainError bellnum(-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]] From c0b96bcdad0c94ddb834031567ceb66c50e20126 Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Fri, 6 Nov 2015 15:10:28 -0500 Subject: [PATCH 2/5] Fix up list of exports --- src/Combinatorics.jl | 2 ++ src/combinations.jl | 2 -- src/permutations.jl | 20 ++++++++++++++++++-- test/runtests.jl | 1 + 4 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/Combinatorics.jl b/src/Combinatorics.jl index be625fa..44fea2b 100644 --- a/src/Combinatorics.jl +++ b/src/Combinatorics.jl @@ -2,6 +2,8 @@ module Combinatorics using Compat, Polynomials, Iterators +import Base: start, next, done, length, eltype + include("numbers.jl") include("factorials.jl") include("combinations.jl") diff --git a/src/combinations.jl b/src/combinations.jl index 7a06e44..87ab390 100644 --- a/src/combinations.jl +++ b/src/combinations.jl @@ -1,8 +1,6 @@ export combinations - #The Combinations iterator -import Base: start, next, done, length immutable Combinations{T} a::T diff --git a/src/permutations.jl b/src/permutations.jl index badabd7..c530fe8 100644 --- a/src/permutations.jl +++ b/src/permutations.jl @@ -1,13 +1,29 @@ +#Permutations + +export + invperm, + ipermute!, + isperm, + levicivita, + nthperm!, + nthperm, + parity, + permutations, + permute!!, + permute! + +#The basic permutaitons iterator + immutable Permutations{T} a::T end +permutations(a) = Permutations(a) + eltype{T}(::Type{Permutations{T}}) = Vector{eltype(T)} length(p::Permutations) = factorial(length(p.a)) -permutations(a) = Permutations(a) - start(p::Permutations) = [1:length(p.a);] function next(p::Permutations, s) perm = [p.a[si] for si in s] diff --git a/test/runtests.jl b/test/runtests.jl index 9cb5de0..94af7fd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,5 +2,6 @@ using Combinatorics using Base.Test include("basic.jl") +include("frombase.jl") include("youngdiagrams.jl") From 9328ef106dd9ad316287053c234e0b4578c81fa2 Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Fri, 6 Nov 2015 15:10:54 -0500 Subject: [PATCH 3/5] Add tests from base --- test/frombase.jl | 67 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 test/frombase.jl diff --git a/test/frombase.jl b/test/frombase.jl new file mode 100644 index 0000000..c612da4 --- /dev/null +++ b/test/frombase.jl @@ -0,0 +1,67 @@ +using Combinatorics +using Base.Test + +import Combinatorics: isperm, invperm, ipermute!, combinations, nthperm, + nthperm!, parity, permute!, permutations, levicivita + +p = shuffle([1:1000;]) +@test isperm(p) +@test all(invperm(invperm(p)) .== p) + +push!(p, 1) +@test !isperm(p) + +a = randcycle(10) +@test ipermute!(permute!([1:10;], a),a) == [1:10;] + +# PR 12785 +let a = 2:-1:1 + @test ipermute!(permute!([1, 2], a), a) == [1, 2] +end + +@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(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 collect(filter(x->(iseven(x[1])),combinations([1,2,3],2))) == Any[[2,3]] + +@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(permutations(0)) == 1 +@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)) + +for n = 0:7, k = 1:factorial(n) + p = nthperm!([1:n;], k) + @test isperm(p) + @test nthperm(p) == k +end + +@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[[]] From 78b31c34fb2f3dd8d45ed1cbde20852de420641d Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Fri, 6 Nov 2015 15:23:24 -0500 Subject: [PATCH 4/5] Reorganize tests Functions from src/foo.jl are now tested in test/foo.jl --- test/basic.jl | 55 ------------------------------------ test/combinations.jl | 17 +++++++++++ test/factorials.jl | 23 +++++++++++++++ test/frombase.jl | 67 -------------------------------------------- test/numbers.jl | 30 ++++++++++++++++++++ test/partitions.jl | 23 +++++++++++++++ test/permutations.jl | 45 +++++++++++++++++++++++++++++ test/runtests.jl | 7 +++-- 8 files changed, 143 insertions(+), 124 deletions(-) delete mode 100644 test/basic.jl create mode 100644 test/combinations.jl create mode 100644 test/factorials.jl delete mode 100644 test/frombase.jl create mode 100644 test/numbers.jl create mode 100644 test/partitions.jl create mode 100644 test/permutations.jl diff --git a/test/basic.jl b/test/basic.jl deleted file mode 100644 index 8be45e1..0000000 --- a/test/basic.jl +++ /dev/null @@ -1,55 +0,0 @@ -using Combinatorics -using Base.Test - -# catalan -@test catalannum(5) == 42 -@test catalannum(30) == parse(BigInt,"3814986502092304") - -#combinations -@test collect(Combinatorics.combinations([])) == [] -@test collect(Combinatorics.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 fibonaccinum(5) == 5 -@test fibonaccinum(101) == parse(BigInt,"573147844013817084101") - -# hyperfactorial -@test hyperfactorial(8) == parse(BigInt,"55696437941726556979200000") - -# 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") - -# 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 bellnum(7) == 877 -@test bellnum(42) == parse(BigInt,"35742549198872617291353508656626642567") -@test_throws DomainError bellnum(-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..b9d0e23 --- /dev/null +++ b/test/factorials.jl @@ -0,0 +1,23 @@ +using Combinatorics +using Base.Test + +# 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/frombase.jl b/test/frombase.jl deleted file mode 100644 index c612da4..0000000 --- a/test/frombase.jl +++ /dev/null @@ -1,67 +0,0 @@ -using Combinatorics -using Base.Test - -import Combinatorics: isperm, invperm, ipermute!, combinations, nthperm, - nthperm!, parity, permute!, permutations, levicivita - -p = shuffle([1:1000;]) -@test isperm(p) -@test all(invperm(invperm(p)) .== p) - -push!(p, 1) -@test !isperm(p) - -a = randcycle(10) -@test ipermute!(permute!([1:10;], a),a) == [1:10;] - -# PR 12785 -let a = 2:-1:1 - @test ipermute!(permute!([1, 2], a), a) == [1, 2] -end - -@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(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 collect(filter(x->(iseven(x[1])),combinations([1,2,3],2))) == Any[[2,3]] - -@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(permutations(0)) == 1 -@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)) - -for n = 0:7, k = 1:factorial(n) - p = nthperm!([1:n;], k) - @test isperm(p) - @test nthperm(p) == k -end - -@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/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..edb0561 --- /dev/null +++ b/test/partitions.jl @@ -0,0 +1,23 @@ +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) diff --git a/test/permutations.jl b/test/permutations.jl new file mode 100644 index 0000000..40434da --- /dev/null +++ b/test/permutations.jl @@ -0,0 +1,45 @@ +using Combinatorics +using Base.Test + +import Combinatorics: invperm, ipermute!, isperm, levicivita, nthperm, nthperm!, parity, permutations, permute! + +p = shuffle([1:1000;]) +@test isperm(p) +@test all(invperm(invperm(p)) .== p) + +push!(p, 1) +@test !isperm(p) + +a = randcycle(10) +@test ipermute!(permute!([1:10;], a),a) == [1:10;] + +# PR 12785 +let a = 2:-1:1 + @test ipermute!(permute!([1, 2], a), a) == [1, 2] +end + +@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 + +for n = 0:7, k = 1:factorial(n) + p = nthperm!([1:n;], k) + @test isperm(p) + @test nthperm(p) == k +end + +@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 94af7fd..7acd6ef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,10 @@ using Combinatorics using Base.Test -include("basic.jl") -include("frombase.jl") +include("numbers.jl") +include("factorials.jl") +include("combinations.jl") +include("permutations.jl") +include("partitions.jl") include("youngdiagrams.jl") From f28bd62852f41ee2ca3e46dbde0c3efe4f733db9 Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Fri, 6 Nov 2015 15:25:11 -0500 Subject: [PATCH 5/5] Migrate functions, docs and tests from Base Also: - clean up readme - expose docs that were buried in comments - get rid of @compat usage --- README.md | 10 +-- REQUIRE | 1 - src/combinations.jl | 11 +-- src/factorials.jl | 30 +++++++-- src/numbers.jl | 16 +++-- src/partitions.jl | 156 ++++++++++++++++++++++++++----------------- src/permutations.jl | 104 +++++++---------------------- src/youngdiagrams.jl | 22 +++--- test/factorials.jl | 9 +++ test/partitions.jl | 5 ++ test/permutations.jl | 22 ++---- 11 files changed, 193 insertions(+), 193 deletions(-) diff --git a/README.md b/README.md index 40ee390..2965b53 100644 --- a/README.md +++ b/README.md @@ -10,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 045ea20..77474e8 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,4 +1,3 @@ julia 0.5- -Compat Polynomials Iterators diff --git a/src/combinations.jl b/src/combinations.jl index 87ab390..ca117cc 100644 --- a/src/combinations.jl +++ b/src/combinations.jl @@ -33,8 +33,8 @@ 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 @@ -43,7 +43,10 @@ function combinations(a, t::Integer) Combinations(a, t) end -#generate combinations of all orders, chaining of order iterators is eager, -#but sequence at each order is lazy + +""" +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 index a4f2307..312835f 100644 --- a/src/factorials.jl +++ b/src/factorials.jl @@ -2,6 +2,7 @@ export derangement, + factorial, subfactorial, doublefactorial, hyperfactorial, @@ -10,7 +11,24 @@ export primorial, multinomial -# The number of permutations of n with no fixed points (subfactorial) +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])) @@ -23,7 +41,7 @@ function doublefactorial(n::Integer) end z = BigInt() ccall((:__gmpz_2fac_ui, :libgmp), Void, - (Ptr{BigInt}, UInt), &z, @compat(UInt(n))) + (Ptr{BigInt}, UInt), &z, UInt(n)) return z end @@ -36,7 +54,7 @@ function multifactorial(n::Integer, m::Integer) end z = BigInt() ccall((:__gmpz_mfac_uiui, :libgmp), Void, - (Ptr{BigInt}, UInt, UInt), &z, @compat(UInt(n)), @compat(UInt(m))) + (Ptr{BigInt}, UInt, UInt), &z, UInt(n), UInt(m)) return z end @@ -46,15 +64,15 @@ function primorial(n::Integer) end z = BigInt() ccall((:__gmpz_primorial_ui, :libgmp), Void, - (Ptr{BigInt}, UInt), &z, @compat(UInt(n))) + (Ptr{BigInt}, UInt), &z, UInt(n)) return z end -#Multinomial coefficient where n = sum(k) +"Multinomial coefficient where n = sum(k)" function multinomial(k...) s = 0 result = 1 - for i in k + @inbounds for i in k s += i result *= binomial(s, i) end diff --git a/src/numbers.jl b/src/numbers.jl index 404e292..354bb83 100644 --- a/src/numbers.jl +++ b/src/numbers.jl @@ -9,7 +9,7 @@ export bellnum, lucasnum, stirlings1 -# Returns the n-th Bell number +"Returns the n-th Bell number" function bellnum(bn::Integer) if bn < 0 throw(DomainError()) @@ -28,7 +28,7 @@ function bellnum(bn::Integer) return list[end] end -# Returns the n-th Catalan number +"Returns the n-th Catalan number" function catalannum(bn::Integer) if bn<0 throw(DomainError()) @@ -44,7 +44,7 @@ function fibonaccinum(n::Integer) end z = BigInt() ccall((:__gmpz_fib_ui, :libgmp), Void, - (Ptr{BigInt}, UInt), &z, @compat(UInt(n))) + (Ptr{BigInt}, UInt), &z, UInt(n)) return z end @@ -56,8 +56,10 @@ function jacobisymbol(a::Integer, b::Integer) (Ptr{BigInt}, Ptr{BigInt}), &ba, &bb) end -#Computes Lassalle's sequence -#OEIS entry A180874 +""" +Computes Lassalle's sequence +OEIS entry A180874 +""" function lassallenum(m::Integer) A = ones(BigInt,m) for n=2:m @@ -79,11 +81,11 @@ function lucasnum(n::Integer) end z = BigInt() ccall((:__gmpz_lucnum_ui, :libgmp), Void, - (Ptr{BigInt}, UInt), &z, @compat(UInt(n))) + (Ptr{BigInt}, UInt), &z, UInt(n)) return z end -# Returns s(n, k), the signed Stirling number of first kind +"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] diff --git a/src/partitions.jl b/src/partitions.jl index 641ca41..8d5b60e 100644 --- a/src/partitions.jl +++ b/src/partitions.jl @@ -1,6 +1,12 @@ #Partitions -export cool_lex, integer_partitions, ncpartitions +export + cool_lex, + integer_partitions, + ncpartitions, + partitions, + prevprod + #nextprod, #integer partitions @@ -11,12 +17,17 @@ end length(p::IntegerPartitions) = npartitions(p.n) -partitions(n::Integer) = IntegerPartitions(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 @@ -77,6 +88,9 @@ 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[] @@ -139,6 +153,9 @@ 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)) @@ -206,6 +223,9 @@ 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) @@ -278,51 +298,59 @@ function nfixedsetpartitions(n::Int,m::Int) return numpart end - -# For a list of integers i1, i2, i3, find the smallest -# i1^n1 * i2^n2 * i3^n3 >= 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 - -# For a list of integers i1, i2, i3, find the largest -# i1^n1 * i2^n2 * i3^n3 <= x -# for integer n1, n2, n3 +#This function is still defined in Base because it is being used by Base.DSP +#""" +#Next integer not less 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 smallest +# i1^n1 * i2^n2 * i3^n3 >= 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")) @@ -362,7 +390,7 @@ function prevprod(a::Vector{Int}, x) end -# Lists the partitions of the number n, the order is consistent with GAP +"Lists the partitions of the number n, the order is consistent with GAP" function integer_partitions(n::Integer) if n < 0 throw(DomainError()) @@ -384,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 index c530fe8..4a54ed5 100644 --- a/src/permutations.jl +++ b/src/permutations.jl @@ -1,23 +1,21 @@ #Permutations export - invperm, - ipermute!, - isperm, levicivita, nthperm!, nthperm, parity, - permutations, - permute!!, - permute! + permutations -#The basic permutaitons iterator +#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)} @@ -46,6 +44,8 @@ function next(p::Permutations, 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")) @@ -66,8 +66,11 @@ function nthperm!(a::AbstractVector, k::Integer) 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) @@ -80,84 +83,21 @@ function nthperm{T<:Integer}(p::AbstractVector{T}) return k end -function invperm(a::AbstractVector) - b = zero(a) # similar vector of zeros - n = length(a) - for i = 1:n - j = a[i] - ((1 <= j <= n) && b[j] == 0) || - throw(ArgumentError("argument is not a permutation")) - b[j] = i - end - b -end - -function isperm(A) - n = length(A) - used = falses(n) - for a in A - (0 < a <= n) && (used[a] $= true) || return false - end - true -end - -function permute!!{T<:Integer}(a, p::AbstractVector{T}) - count = 0 - start = 0 - while count < length(a) - ptr = start = findnext(p, start+1) - temp = a[start] - next = p[start] - count += 1 - while next != start - a[ptr] = a[next] - p[ptr] = 0 - ptr = next - next = p[next] - count += 1 - end - a[ptr] = temp - p[ptr] = 0 - end - a -end - -permute!(a, p::AbstractVector) = permute!!(a, copy!(similar(p), p)) - -function ipermute!!{T<:Integer}(a, p::AbstractVector{T}) - count = 0 - start = 0 - while count < length(a) - start = findnext(p, start+1) - temp = a[start] - next = p[start] - count += 1 - while next != start - temp_next = a[next] - a[next] = temp - temp = temp_next - ptr = p[next] - p[next] = 0 - next = ptr - count += 1 - end - a[next] = temp - p[next] = 0 - end - a -end - -ipermute!(a, p::AbstractVector) = ipermute!!(a, copy!(similar(p), p)) +# 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. -# 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. -# Returns 1 is the permutarion is even, -1 if it is odd and 0 otherwise. +""" +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) @@ -187,8 +127,10 @@ function levicivita{T<:Integer}(p::AbstractVector{T}) 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. +""" +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")) 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/factorials.jl b/test/factorials.jl index b9d0e23..a9342bb 100644 --- a/test/factorials.jl +++ b/test/factorials.jl @@ -1,6 +1,15 @@ 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") diff --git a/test/partitions.jl b/test/partitions.jl index edb0561..7f90ddf 100644 --- a/test/partitions.jl +++ b/test/partitions.jl @@ -21,3 +21,8 @@ using Base.Test # 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 index 40434da..e4275b1 100644 --- a/test/permutations.jl +++ b/test/permutations.jl @@ -1,22 +1,7 @@ using Combinatorics using Base.Test -import Combinatorics: invperm, ipermute!, isperm, levicivita, nthperm, nthperm!, parity, permutations, permute! - -p = shuffle([1:1000;]) -@test isperm(p) -@test all(invperm(invperm(p)) .== p) - -push!(p, 1) -@test !isperm(p) - -a = randcycle(10) -@test ipermute!(permute!([1:10;], a),a) == [1:10;] - -# PR 12785 -let a = 2:-1:1 - @test ipermute!(permute!([1, 2], a), a) == [1, 2] -end +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']] @@ -26,11 +11,16 @@ end @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])