Combinatorics

catalan(n::fmpz; alg="binomial")
catalan(n::Integer; alg="binomial")

The n-th Catalan number. This counts a gazillion of things, see OEIS for more information. There are two algorithms implemented:

1. "binomial" (default): uses a simple formula with binomial coefficients.
2. "iterative": uses an iterative computation.

The binomial computation is much faster:

julia> @time x=catalan( ZZ(10)^5 , alg="binomial");
 0.007727 seconds (9 allocations: 95.750 KiB)

julia> @time x=catalan( ZZ(10)^5 , alg="iterative");
 1.572488 seconds (1.01 M allocations: 2.333 GiB, 1.36% gc time)

num_partitions(n::fmpz)
num_partitions(n::Integer)

The number of integer partitions of the integer $n ≥ 0$. Uses the function from FLINT, which is really fast.

For more information on these numbers, see OEIS.

num_partitions(n::fmpz, k::fmpz)
num_partitions(n::Integer, k::Integer)

The number of integer partitions of the integer $n ≥ 0$ into $k ≥ 0$ parts. The implementation uses a recurrence relation.

For more information on these numbers, see OEIS.

lucas(n::fmpz)
lucas(n::Integer)

The n-th Lucas number. For more information on these numbers, see OEIS. The implementation is a wrapper to the function in GMP.

stirling1(n::fmpz, k::fmpz)
stirling1(n::Integer, k::Integer)

The Stirling number $S_1(n,k)$ of the first kind. The absolute value of $S_1(n,k)$ counts the number of permutations of $n$ elements with $k$ disjoint cycles. The implementation is a wrapper to the function in FLINT.

For more information on these numbers, see OEIS.

stirling2(n::fmpz, k::fmpz)
stirling2(n::Integer, k::Integer)

The Stirling number $S_2(n,k)$ of the second kind. This counts the number of partitions of an $n$-set into $k$ non-empty subsets. The implementation is a wrapper to the function in FLINT.

For more information on these numbers, see OEIS.

Partition{T} <: AbstractArray{T,1}

A partition of an integer n ≥ 0 is a decreasing sequence λ=(λ₁,…,λᵣ) of positive integers λᵢ whose sum is equal to n. The λᵢ are called the parts of the partition. We encode a partition as an array with elements λᵢ. To be able to conceptually work with partitions we have implemented an own type Partition{T} as subtype of AbstractArray{T,1}. All functions for arrays then also work for partitions. You may increase performance by using smaller integer types, see the example below. For efficiency, the Partition constructor does not check whether the given array is in fact a partition, i.e. a decreasing sequence.

For more general information on partitions, check out Wikipedia.

Example

julia> P=Partition([3,2,1]) #The partition 3+2+1 of 6
julia> sum(P) #The sum of the parts.
6
julia> P[1] #First component
3
julia> P=Partition(Int8[3,2,1]) #Same partition but using 8 bit integers

Remarks

• Usually, |λ| ≔ n is called the size of λ. In Julia, the function size for arrays already exists and returns the dimension of an array. Instead, you can use sum to get the sum of the parts.

• There is no performance impact by using an own type for partitions rather than simply using arrays—I've tested this. Julia is great. The implementation of a subtype of AbstractArray is explained in the Julia documentation.

partitions(n::Integer)

A list of all partitions of an integer n ≥ 0, produced in lexicographically descending order. This ordering is like in SAGE, but opposite to GAP. You can apply reverse() to reverse the order. As usual, you may increase performance by using smaller integer types.

The algorithm used is the algorithm ZS1 by A. Zoghbi and I. Stojmenovic, "Fast algorithms for generating integer partitions", Int. J. Comput. Math. 70 (1998), no. 2, 319–332.

Example

julia> partitions(Int8(90)) #Using 8-bit integers

partitions(m::Integer, n::Integer, l1::Integer, l2::Integer; z=0)

A list of all partitions of an integer m ≥ 0 into n ≥ 0 parts with lower bound l1 ≥ 0 and upper bound l2 ≥ l1 for the parts. There are two choices for the parameter z:

• z=0: no further restriction (default);
• z=1: only distinct parts.

The partitions are produced in decreasing order.

The algorithm used is "parta" by W. Riha and K. R. James, "Algorithm 29. Efficient Algorithms for Doubly and Multiply Restricted Partitions" (1976). De-gotoed from ALGOL code by Elisa!

partitions(m::Integer, n::Integer)

All partitions of an integer m ≥ 0 into n ≥ 1 parts (no further restrictions). This simply calls partitions(m,n,1,m,z=0).

partitions(mu::Array{Integer,1}, m::Integer, v::Array{Integer,1}, n::Integer)

All partitions of an integer m >= 0 into n >= 1 parts, where each part is an element in v and each v[i] occurs a maximum of mu[i] times. The partitions are produced in decreasing order.

The algorithm used is a de-gotoed version of "partb" by W. Riha and K. R. James, "Algorithm 29. Efficient Algorithms for Doubly and Multiply Restricted Partitions" (1976). Note: The original Algorithm lead to BoundsErrors, since r could get smaller than 1. Furthermore x and y are handled as Arrays with an infinite length. After finding all valid partitions, the algorithm will continue searching for partitions of length n+1. We thus had to add a few additional checks and interruptions.

dominates(λ::Partition, μ::Partition)

The dominance order on partitions is the partial order ⊵ defined by λ ⊵ μ if and only if λ₁ + … + λᵢ ≥ μ₁ + … + μᵢ for all i. This function returns true if λ ⊵ μ.

For more information see Wikipedia.

conjugate(λ::Partition{T}) where T<:Integer

The conjugate of a partition is obtained by considering its Young diagram (see Tableau) and then flipping it along its main diagonal.

For more information see Wikipedia.

getelement(P::Partition, i::Int)

Sometimes in algorithms for partitions it is convenient to be able to access parts beyond the length of the partition, and then you want to get zero instead of an error. This function is a shortcut for

return (i>length(P.p) ? 0 : getindex(P.p,i))

If you are sure that P[i] exists, use getindex because this will be faster.

ascending_partitions(n::Integer;alg="ks")

Instead of encoding a partition of an integer n ≥ 0 as a descending sequence (which is our convention), one can also encode it as an ascending sequence. In the papers below it is claimed that generating the list of all ascending partitions is more efficient than generating descending ones. To test this, I have implemented the algorithms:

1. "ks" (default) is the algorithm AccelAsc (Algorithm 4.1) by J. Kelleher and B. O'Sullivan, "Generating All Partitions: A Comparison Of Two Encodings", https://arxiv.org/pdf/0909.2331.pdf, May 2014.
2. "m" is Algorithm 6 by M. Merca, "Fast Algorithm for Generating Ascending Compositions", J. Math Model. Algor. (2012) 11:89–104. This is similar to "ks".

The ascending partitions are given here as arrays, not of type Partition since these are descending by our convention. I am using "ks" as default since it looks slicker and I believe there is a tiny mistake in the publication of "m" (which I fixed).

Comparison

I don't see a significant speed difference to the descending encoding:

julia> @btime partitions(Int8(90));
  3.376 s (56634200 allocations: 6.24 GiB)

julia> @btime ascending_partitions(Int8(90),alg="ks");
  3.395 s (56634200 allocations: 6.24 GiB)

julia> @btime ascending_partitions(Int8(90),alg="m");
  3.451 s (56634200 allocations: 6.24 GiB)

Multipartition{T} <: AbstractArray{Partition{T},1}

Multipartitions are generalizations of partitions. An r-component multipartition of an integer n is an r-tuple of partitions λ¹, λ², …, λʳ where each λⁱ is a partition of some integer nᵢ ≥ 0 and the nᵢ sum to n. As for partitions, we have implemented an own type Multipartition{T} which is a subtype of AbstractArray{Partition{T},1}. As with partitions, you can can use smaller integer types to increase performance.

Example

julia> P=Multipartition( [[2,1], [], [3,2,1]] )
julia> sum(P)
9
julia> P[2]
Int64[]
julia> Multipartition( Array{Int8,1}[[2,1], [], [3,2,1]] ) #Using 8-bit integers

multipartitions(n::T, r::Integer) where T<:Integer

A list of all r-component multipartitions of n.

Example

julia> multipartitions(Int8(3),2) #Using 8-bit integers

sum(P::Multipartition{T}) where T<:Integer

If P is a multipartition of the integer n, this function returns n.

Multiset_partition{T}

Multiset-partitions are generalizations of partitions. An r-component multiset-partition of an integer $n$ is a multiset(a set where each element can be contained multiple times) of partitions $λ¹, λ², …, λʳ$ where each $λⁱ$ is a partition of some integer $nᵢ ≥ 1$ and the $nᵢ$ sum to $n$. As for partitions, we have implemented an own type Multiset_partition{T}. As with partitions, you can can use smaller integer types to increase performance.

Example

julia> P = Multiset_partition( [2,1], [4], [3,2,1] )
{[2, 1], [4], [3, 2, 1]}
julia> sum(P)
13
julia> Multiset_partition( Array{Int8,1}[[2,1], [4], [3,2,1]] ) #Using 8-bit integers
{Int8[2, 1], Int8[4], Int8[3, 2, 1]}

Since Multiset-partitions are unordered sets, you can't call an explicit element, however, you can iterate over a Multiset_partition.

Example

julia> MSP = Multiset_partition( [2,1], [4], [3,2,1] )
{[2, 1], [4], [3, 2, 1]}
julia> for p in MSP println(p) end
[2, 1]
[4]
[3, 2, 1]

multiset_partitions(n::T) where T<:Integer

A list of all multiset_partitions of an integer $n ⋝ 0$.

The performance will suffer by casting $n$ into a smaller integer type, e.g.

multiset_partitions(Int8(20))

multiset_partitions(p::Partition{T})  where T<:Integer

A list of all possible multiset_partitions of a Partition, by regrouping its parts into Partitions.

The algorithm used is the algorithm M by , "The Art of Computer Programming - Volume 4A, Combinatorial Algorithms, Part 1" by Donald E. Knuth(2011), 429–430. De-gotoed, index-shifted and generalized.

multiset_partitions(n::T, r::Integer) where T<:Integer

A list of all multiset_partitions of an integer $n ⋝ 0$ into $r ⋝ 1$ parts.

multiset_partitions(p::Partition{T}, r::Integer) where T<:Integer

A list of all possible $r$-restricted multiset_partitions of a Partition, by regrouping its parts into Partitions.

The algorithm used is a version of the algorithm M by , "The Art of Computer Programming - Volume 4A, Combinatorial Algorithms, Part 1" by Donald E. Knuth, 429–430 http://www.cs.utsa.edu/~wagner/knuth/fasc3b.pdf. De-gotoed, index-shifted and generalized.

Tableau{T} <: AbstractArray{AbstractArray{T,1},1}

A Young diagram is a diagram of finitely many empty "boxes" arranged in left-justified rows, with the row lengths in non-increasing order. The box in row i and and column j has the coordinates (i,j). Listing the number of boxes in each row gives a partition λ of a non-negative integer n (the total number of boxes of the diagram). The diagram is then said to be of shape λ. Conversely, one can associate to any partition λ a Young diagram in the obvious way, so Young diagrams are just another way to look at partitions.

A Young tableau of shape λ is a filling of the boxes of the Young diagram of λ with elements from some set. After relabeling we can (and will) assume that we fill from a set of integers from 1 up to some number, which in applications is often equal to n. We encode a tableau as an array of arrays and we have implemented an own type Tableau{T} as subtype of AbstractArray{AbstractArray{T,1},1} to work with tableaux. As for partitions, you may increase performance by casting into smaller integer types, e.g.

For efficiency, we do not check whether the given array is really a tableau, i.e. whether the structure of the array defines a partition.

For more information see Wikipedia.

Example

Tab=Tableau([[1,2,3],[4,5],[6]])
Tab=Tableau(Array{Int8,1}[[2,1], [], [3,2,1]]) #Using 8 bit integers

shape(Tab::Tableau{T})

Returns the shape of a tableau.

weight(Tab::Tableau)

The weight of a tableau is the number of times each number appears in the tableau. The return value is an array whose i-th element gives the number of times the integer i appears in the tableau.

reading_word(Tab::Tableau)

The reading word of a tableau is the word obtained by concatenating the fillings of the rows, starting from the bottom row. The word is here returned as an array.

Example

julia> reading_word(Tableau([ [1,2,3] , [4,5] , [6] ]))
6-element Array{Int64,1}:
 6
 4
 5
 1
 2
 3

is_semistandard(Tab::Tableau)

A tableau is called semistandard if the entries weakly increase along each row and strictly increase down each column.

semistandard_tableaux(shape::Partition{T}, max_val=sum(shape)::Integer) where T<:Integer

Returns a list of all semistandard tableaux of given shape and filling elements bounded by max_val. By default, max_val is equal to the sum of the shape partition (the number of boxes in the Young diagram). The list of tableaux is in lexicographic order from left to right and top to bottom.

semistandard_tableaux(shape::Partition{T}, max_val=sum(shape)::Integer) where T<:Integer

Shortcut for semistandard_tableaux(Partition(shape),max_val).

semistandard_tableaux(box_num::T, max_val=box_num::T) where T<:Integer

Returns a list of all semistandard tableaux consisting of box_num boxes and filling elements bounded by max_val.

semistandard_tableaux(s::Array{T,1}, weight::Array{T,1}) where T<:Integer

Returns a list of all semistandard tableaux with shape s and given weight. This requires that sum(s) = sum(weight).

semistandard_tableaux(s::Partition{T}, weight::Partition{T}) where T<:Integer

Same as for semistandard_tableaux(s::Array{T,1}, weight::Array{T,1}) where T<:Integer.

is_standard(Tab::Tableau)

A tableau is called standard if it is semistandard and the entries are in bijection with 1,…n, where n is the number of boxes.

standard_tableaux(s::Partition)

Returns a list of all standard tableaux of a given shape.

standard_tableaux(s::Array{Integer,1})

Shortcut for standard_tableaux(Partition(s)).

standard_tableaux(n::Integer)

Returns a list of all standard tableaux with n boxes.

hook_length(lambda::Partition, i::Integer, j::Integer)

Consider the Young diagram of a partition λ. The hook length of a box, is the number of boxes to the right in the same row + the number of boxes below in the same column + 1. The function returns the hook length of the box with coordinates (i,j). The functions assumes that the box exists.

hook_length(Tab::Tableau, i::Integer, j::Integer)

Shortcut for hook_length(shape(Tab),i,j).

hook_lengths(lambda::Partition)

Returns the tableau of shape λ in which the entry at position (i,j) is equal to the hook length of the corresponding box.

num_standard_tableaux(lambda::Partition)

Returns the number $f^\lambda$ of standard tableaux of shape λ using the hook length formula

\[f^{\lambda} = \frac{n!}{\prod_{i,j} h_\lambda(i,j)} \;,\]

where the product is taken over all boxes in the Young diagram of $\lambda$ and $h_\lambda$ denotes the hook length of the box (i,j).

For my information, see Wikipedia.

schensted(sigma::Array{Integer,1})

The Robinson–Schensted correspondence is a bijection between permutations and pairs of standard Young tableaux of the same shape. For a permutation sigma (given as an array), this function performs the Schnested algorithm and returns the corresponding pair of standard tableaux (the insertion and recording tableaux).

For more information, see Wikipedia.

Example

julia> schensted([3,1,6,2,5,4])
(AbstractArray{Int64,1}[[1, 2, 4], [3, 5], [6]], AbstractArray{Int64,1}[[1, 3, 5], [2, 4], [6]])

bump!(Tab::Tableau, x::Int)

Inserts the integer x into the tableau Tab according to the bumping algorithm by applying the Schensted insertion.

bump!(Tab::Tableau, x::Integer, Q::Tableau, y::Integer)

Inserts x into Tab according to the bumping algorithm by applying the Schensted insertion. Traces the change with Q by inserting y at the same Position in Q as x in Tab.

kostka_polynomial(λ::Partition{T}, μ::Partition{T})
kostka_polynomial(λ::Array{Int,1}, μ::Array{Int,1})

The (one-variable) Kostka polymial $K_{λμ}(t)$ associated to partitions λ and μ can be defined as

\[K_{λμ}(t) = ∑_{T ∈ SSYT(λ,μ)} t^{charge(T)} ∈ ℕ[t] \;,\]

where SSYT(λ,μ) is the set of all semistandard Young tableaux of shape λ and weight μ, and charge(T) is the charge of the tableau T. The Kostka polynomials relate the Hall–Littlewood polynomials $P_μ$ to the Schur polynomials $s_λ$ via

\[s_λ(x_1,…,x_n) = \sum_μ K_{λμ}(t)P_μ(x_1,…,x_n;t)\]

This function returns the Kostka polynomial $K_{λμ}(t)$ as an fmpz_poly over ZZ in t.

Example

julia> kostka_polynomial([4,2,1],[3,2,1,1])
t^3 + 2*t^2 + t

Algorithm

The computation here is not based on the above formula but on an explicit description due to Kirillov–Reshetikhin in "The Bethe ansatz and the combinatorics of Young tableaux", J. Sov. Math. 41 (1988) 925., which is summarized by Dorey–Tonga–Turner in "A matrix model for WZW" (Apendix B). Namely:

\[K_{λ,μ}(t)=∑_{\{v\}}∏_{K=1}^{l(λ)-1}∏_{n≥1} \begin{bmatrix} ℙ_n^{(K)}+m_n(v^{(K)})\\ m_n(v^{(K)}) \end{bmatrix}_t \;,\]

where the sum is over all admissible configurations $\{v\}$ i.e. sequences of partitions $v^{(K)}$ with

\[v^{(0)}=μ \text{ ,\hspace{2mm} } |v^{(K)}|=∑_{j≥K+1}λ_j \text{\hspace{5mm} and \hspace{5mm}} ℙ_n^{(K)}≥0 \text{\hspace{2mm} for all \hspace{1mm}} n>0, K=0,1,…l(λ)\]

\[\begin{aligned} & ℙ_n^{(K)} := ∑_{j≥1}\left[\min\left(n,v_j^{(K+1)}\right) - 2⋅\min\left(n,v_j^{(K)}\right) + \min\left(n,v_j^{(K-1)}\right)\right] \\ & c(v) := ∑_{i≥1}(i-1)μ_i + ∑_{K=1}^{l(λ)-1}\left(𝕄\left[v^{(K)}, v^{(K)}\right] - 𝕄\left[v^{(K)}, v^{(K-1)}\right]\right) \\ & 𝕄[ρ,κ] := ∑_{i,j≥1} \min(ρ_i,κ_j) \end{aligned}\]

Here, $\left[\genfrac{[}{]}{0pt}{0}{m}{n} \right]_t$ is the Gaussian binomial coefficient.

charge(config::Array{Partition{T},1})

The charge $c$ of an admissible configuration $v:=$ config, is defined by:

\[\begin{aligned} & c(v) := ∑_{i≥1}(i-1)v^{(0)}_i + ∑_{K=1}^{l(λ)-1}\left(𝕄\left[v^{(K)}, v^{(K)}\right] - 𝕄\left[v^{(K)}, v^{(K-1)}\right]\right) \\ & 𝕄[ρ,κ] := ∑_{i,j≥1} \min(ρ_i,κ_j) \end{aligned}\]

charge(Tab::Tableau)

returns the charge of a Tableau Tab which is defined by the charge of it´s reading word.

charge(Tab):=charge(reading_word(Tab))

charge(word::Array{Int,1},standard=false::Bool)

This returns the charge of the Tableau corresponding to the reading word word.

We call a word standard, if each of its letters or numbers only appears once in it. If you are shure that your word is standard, you can set this argument to $true$, to improve the efficiency.

This Algorithm is based on the Algorithm following Example 7.3 in "Hall-Littlewood Functions and Kostka-Foulkes Polynomials in Representation Theory", J. Désarménien, B. Leclerc and J.-Y. Thibon.

schur_polynomial(λ::Partition{T}, n=sum(λ)::Int) where T<:Integer
schur_polynomial(λ::Partition{T}, R::FmpzMPolyRing, n=sum(λ)::Int) where T<:Integer
schur_polynomial(λ::Partition{T}, x::Array{fmpz_mpoly,1}) where T<:Integer

Returns the Schur polynomial $s_λ(x₁,x₂,...,xₙ)$ in n variables, as a Multivariate Polynomial.

If neither R nor x are given, the Schur polynomial will be over PolynomialRing(ZZ,["x1","x2",...,"xn"]).

Example

julia> R,x = PolynomialRing(ZZ, ["a","b","c"])
(Multivariate Polynomial Ring in a, b, c over Integer Ring, fmpz_mpoly[a, b, c])
julia> schur_polynomial(Partition([2,1]),[x[1],x[2]])
a^2*b + a*b^2
julia> schur_polynomial(Partition([2,1]),R)
a^2*b + a^2*c + a*b^2 + 2*a*b*c + a*c^2 + b^2*c + b*c^2
julia> schur_polynomial(Partition([2]))
x1^2 + x1*x2 + x2^2

Algorithm

We use two different Algorithms, depending on the size of the input. The Combinatorial Algorithm is used for Partitions of small Integers, or if $n ≥ 10$. In the other cases we use Cauchy's bialternant formula.

Combinatorial Algorithm

\[s_λ:=∑_T x₁^{m₁}…xₙ^{mₙ}\]

where the sum is taken over all semistandard tableaux $T$ of shape $λ$, and $mᵢ$ gives the weight of $i$ in $T$.

Cauchy's bialternant formula

\[s_λ(x₁,…,xₙ) = ∏_{1 ≤ i < j ≤ n} (x_i-x_j)^{-1} ⋅ \begin{vmatrix} x_1^{λ₁+n-1} & x_2^{λ_1+n-1} & … & x_n^{λ_1+n-1} \\ x_1^{λ_2+n-2} & x_2^{λ_2+n-2} & … & x_n^{λ_2+n-2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ x_1^{λ_n} & x_2^{λ_n} & … & x_n^{λ_n} \end{vmatrix}\]

quantum_number(n::Int, q::RingElem)

For an integer n ≥ 0 and an invertible element q of a ring R, the quantum integer $[n]_q \in R$ is for n ≥ 0 defined as $[n]_q = \sum_{i=0}^{n-1} q^{n-(2i+1)}$ and for n < 0 as $[n]_q = -[-n]_q$.

quantum_number(n::Int)

The quantum number $[n]_q$ where q is the indeterminate of the Laurent polynomial ring $\mathbb{Z}[q,q^{-1}]$ in one variable over the integers.

quantum(n::Int, q::RingElem)

Well This is a shortcut for quantum_number(n,q).

quantum(n::Int)

This is a shortcut for quantum_number(n).

gaussian_binomial(m::Int, r::Int, q::RingElem)

The Gaussian binomial coefficient $\begin{bmatrix} m \\ r \end{bmatrix}_q$ , also known as q-binomial coefficient is defined by:

\[\begin{bmatrix} m \\ r \end{bmatrix}_q = \begin{cases}{\frac{(1-q^m)(1-q^{m-1})⋯(1-q^{m-r+1})}{(1-q)(1-q^2)⋯(1-q^r)}} & r≤m \\ 0 & r>m \end{cases} \;\]