Combinations

To make a matrix with all combinations of the elements of \(S\) (of size \(n\)) taken \(r\) at a time:

S <- letters[1:5]
n <- length(S)
r <- 3
result <- t(combn(S,r))
result
##       [,1] [,2] [,3]
##  [1,] "a"  "b"  "c" 
##  [2,] "a"  "b"  "d" 
##  [3,] "a"  "b"  "e" 
##  [4,] "a"  "c"  "d" 
##  [5,] "a"  "c"  "e" 
##  [6,] "a"  "d"  "e" 
##  [7,] "b"  "c"  "d" 
##  [8,] "b"  "c"  "e" 
##  [9,] "b"  "d"  "e" 
## [10,] "c"  "d"  "e"
apply(result,1,function(x) paste0(x,collapse="")) # just return as a vector of strings
##  [1] "abc" "abd" "abe" "acd" "ace" "ade" "bcd" "bce" "bde" "cde"

The number of possible combinations is \[C(n,r) = \frac{n!}{r!(n-r)!}\]

This number is computed in R by the choose function:

choose(5,0)
## [1] 1
choose(5,1)
## [1] 5
choose(5,2)
## [1] 10
choose(5,3)
## [1] 10
choose(5,4)
## [1] 5
choose(5,5)
## [1] 1

Some recursive relations:

Combinations, among many other things, give the coefficients for the binomial expression:

\[(a+b)^n = \sum_{i=0}^n C(n,i) a^{n-i}b^i\]

and thus combinations give us Pascal Triangle:

pascal <- matrix(rep(NA,100),ncol=10)
for(i in 0:10)
  for(j in 0:i)
    pascal[i,j] <- choose(i,j)
print(pascal, na.print=" " )
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1                                              
##  [2,]    2    1                                         
##  [3,]    3    3    1                                    
##  [4,]    4    6    4    1                               
##  [5,]    5   10   10    5    1                          
##  [6,]    6   15   20   15    6    1                     
##  [7,]    7   21   35   35   21    7    1                
##  [8,]    8   28   56   70   56   28    8    1           
##  [9,]    9   36   84  126  126   84   36    9    1      
## [10,]   10   45  120  210  252  210  120   45   10     1

Permutations

To make the same but where the order is relevant, ie, a permutation (if \(m\) is equal to the length of \(S\) we get all possible arrangements):

library(gtools)

result <- permutations(n, r, S)
apply(result,1,function(x) paste0(x,collapse=""))
##  [1] "abc" "abd" "abe" "acb" "acd" "ace" "adb" "adc" "ade" "aeb" "aec"
## [12] "aed" "bac" "bad" "bae" "bca" "bcd" "bce" "bda" "bdc" "bde" "bea"
## [23] "bec" "bed" "cab" "cad" "cae" "cba" "cbd" "cbe" "cda" "cdb" "cde"
## [34] "cea" "ceb" "ced" "dab" "dac" "dae" "dba" "dbc" "dbe" "dca" "dcb"
## [45] "dce" "dea" "deb" "dec" "eab" "eac" "ead" "eba" "ebc" "ebd" "eca"
## [56] "ecb" "ecd" "eda" "edb" "edc"

The number of possible permutations is \[P(n,r) = \frac{n!}{(n-r)!}\]

Some recursive relations:

A permutation using all the elements of \(S\) can be seen as a bijective mapping \(S \rightarrow S\). There are \(P(n,n)=n!\) such mappings. In this context a fixed point is when an element is mapped into itself. Mappings with no fixed points are called derangements (see below).

No restrictions

To return all possibilities of \(r\) values without any restrictions from the \(n\) elements of \(S\) (an element can be repeated):

result <- as.matrix(expand.grid(lapply(numeric(r), function(x) S)), ncol=r)
apply(result,1,function(x) paste0(x,collapse="")) 
##   [1] "aaa" "baa" "caa" "daa" "eaa" "aba" "bba" "cba" "dba" "eba" "aca"
##  [12] "bca" "cca" "dca" "eca" "ada" "bda" "cda" "dda" "eda" "aea" "bea"
##  [23] "cea" "dea" "eea" "aab" "bab" "cab" "dab" "eab" "abb" "bbb" "cbb"
##  [34] "dbb" "ebb" "acb" "bcb" "ccb" "dcb" "ecb" "adb" "bdb" "cdb" "ddb"
##  [45] "edb" "aeb" "beb" "ceb" "deb" "eeb" "aac" "bac" "cac" "dac" "eac"
##  [56] "abc" "bbc" "cbc" "dbc" "ebc" "acc" "bcc" "ccc" "dcc" "ecc" "adc"
##  [67] "bdc" "cdc" "ddc" "edc" "aec" "bec" "cec" "dec" "eec" "aad" "bad"
##  [78] "cad" "dad" "ead" "abd" "bbd" "cbd" "dbd" "ebd" "acd" "bcd" "ccd"
##  [89] "dcd" "ecd" "add" "bdd" "cdd" "ddd" "edd" "aed" "bed" "ced" "ded"
## [100] "eed" "aae" "bae" "cae" "dae" "eae" "abe" "bbe" "cbe" "dbe" "ebe"
## [111] "ace" "bce" "cce" "dce" "ece" "ade" "bde" "cde" "dde" "ede" "aee"
## [122] "bee" "cee" "dee" "eee"

This is a vector with \(r^n\) elements.

The next table shows how to count the total possibilities for the following conditions: \[ \begin{array} {|c|c|c|} \hline \text{Order Counts?} & \text{Repetition?} & \text{Expression} \\ \hline Yes & Yes & n^r \\ Yes & No & P(n,r) \\ No & Yes & C(n+r-1,r) \\ No & No & C(n,r) \\ \hline \end{array} \]

Circular Permutations

A circular permutation is a circular arrangement of elements for which the order of the elements must be taken into account.

The number of circular permutations of \(r\) elements taken from a set of \(n\) elements is:

\[\frac{P(n,r)}{r}\]

with \(2 \leq r \leq n\).

To produce these circular permutations:

library(binhf) # uses: shift()

result <- permutations(n, r, S)

# receive a vector of comparables, shift the vector so that
# the minimum is at the head of the vector
shift.to.minimum <- function(row) {
  shift(row, which(row==min(row))-1, dir='left')
}

result <- unique(t(apply(result, 1, shift.to.minimum)))
head(result,12)
##       [,1] [,2] [,3]
##  [1,] "a"  "b"  "c" 
##  [2,] "a"  "b"  "d" 
##  [3,] "a"  "b"  "e" 
##  [4,] "a"  "c"  "b" 
##  [5,] "a"  "c"  "d" 
##  [6,] "a"  "c"  "e" 
##  [7,] "a"  "d"  "b" 
##  [8,] "a"  "d"  "c" 
##  [9,] "a"  "d"  "e" 
## [10,] "a"  "e"  "b" 
## [11,] "a"  "e"  "c" 
## [12,] "a"  "e"  "d"

Necklace Permutations

These are 3D circular permutations, ie, they can be flipped. So, there are half of them

\[\frac{P(n,r)}{2r}\]

with \(3 \leq r \leq n\).

result <- permutations(n, r, S)

# shift the vector and its reverse, and select the smallest one (using lexicographic order)
# pre: elements must be pastable into a string
flip.and.shift.to.minimum <- function(row) {
  row1 <- shift.to.minimum(row)
  row2 <- shift.to.minimum(rev(row))
  if (paste(row1,collapse="") < paste(row2,collapse=""))
    row1
  else
    row2
}

result <- unique(t(apply(result, 1, flip.and.shift.to.minimum)))
result
##       [,1] [,2] [,3]
##  [1,] "a"  "b"  "c" 
##  [2,] "a"  "b"  "d" 
##  [3,] "a"  "b"  "e" 
##  [4,] "a"  "c"  "d" 
##  [5,] "a"  "c"  "e" 
##  [6,] "a"  "d"  "e" 
##  [7,] "b"  "c"  "d" 
##  [8,] "b"  "c"  "e" 
##  [9,] "b"  "d"  "e" 
## [10,] "c"  "d"  "e"

Multisets

A multiset (or bag) is a generalization of the notion of a set in which members are allowed to appear more than once.

If we want to arrange multisets we can do something like this:

S <- c(1,2,3)
n <- c(2,3,1)  # we want two 1s, three 2s, and one 3

S1 <- rep(S,n)
S1
## [1] 1 1 2 2 2 3

result <- unique(permutations(length(S1), length(S1), S1, set=FALSE))
head(result,20)
##       [,1] [,2] [,3] [,4] [,5] [,6]
##  [1,]    1    1    2    2    2    3
##  [2,]    1    1    2    2    3    2
##  [3,]    1    1    2    3    2    2
##  [4,]    1    1    3    2    2    2
##  [5,]    1    2    1    2    2    3
##  [6,]    1    2    1    2    3    2
##  [7,]    1    2    1    3    2    2
##  [8,]    1    2    2    1    2    3
##  [9,]    1    2    2    1    3    2
## [10,]    1    2    2    2    1    3
## [11,]    1    2    2    2    3    1
## [12,]    1    2    2    3    1    2
## [13,]    1    2    2    3    2    1
## [14,]    1    2    3    1    2    2
## [15,]    1    2    3    2    1    2
## [16,]    1    2    3    2    2    1
## [17,]    1    3    1    2    2    2
## [18,]    1    3    2    1    2    2
## [19,]    1    3    2    2    1    2
## [20,]    1    3    2    2    2    1

The number of multisets of order \(\{n_1,n_2,\ldots,n_k\}\) is given by

\[\frac{n!}{n_1!n_2!\ldots n_k!}\]

where \(n_i\) means how many i-th elements are there, and \(n = \sum_i n_i\).

This expression can also be used to calculate the coefficients of the multinomial expressions like

\[(a_1 + a_2 + \ldots + a_m)^n = \sum_{k_1+k_2+\ldots+k_m=n} \frac{n!}{k_1!k_2!\ldots k_k!}\]

Eg: The coefficient of \(ac^2\) from \((a+b+c)^3\) is

\[(a+b+c)^3 = \ldots + \frac{3!}{1!0!2!} ac^2 + \ldots\]

Derangements

A derangement is a permutation of the elements of a set such that none of the elements appear in their original position.

# number of derangements, also called subfactorial !n
# pre: n>=0
dn <- function(n) {
  ifelse (n<2, 1-n, round(factorial(n)/exp(1)))
}

print(data.frame(n=0:12, dn=dn(0:12)), row.names=FALSE, quote=FALSE)
##   n        dn
##   0         1
##   1         0
##   2         1
##   3         2
##   4         9
##   5        44
##   6       265
##   7      1854
##   8     14833
##   9    133496
##  10   1334961
##  11  14684570
##  12 176214841

This following relation holds between the subfactorial and factorial:

\[\lim_{n \rightarrow \infty} \frac{!n}{n!} = \frac{1}{e}\]

which is the probability that a randomly selected permutation is a derangement.

options(digits=9)
xs <- 1:14
print(data.frame(n=xs, dn=dn(xs), fact=factorial(xs), ratio=dn(xs)/factorial(xs)), row.names=FALSE, quote=FALSE)
##   n          dn        fact       ratio
##   1           0           1 0.000000000
##   2           1           2 0.500000000
##   3           2           6 0.333333333
##   4           9          24 0.375000000
##   5          44         120 0.366666667
##   6         265         720 0.368055556
##   7        1854        5040 0.367857143
##   8       14833       40320 0.367881944
##   9      133496      362880 0.367879189
##  10     1334961     3628800 0.367879464
##  11    14684570    39916800 0.367879439
##  12   176214841   479001600 0.367879441
##  13  2290792932  6227020800 0.367879441
##  14 32071101049 87178291200 0.367879441
1/exp(1)
## [1] 0.367879441
options(digits=7)

The next script computes all possible derangements:

# computing all derangements, ie, 
library(gtools)
r <- 5
S <- 1:r

results <- permutations(r, r, S)                             # get all permutations
results <- results[apply(results,1,function(x) !any(x==S)),] # filter those with fixed points
head(results,12)
##       [,1] [,2] [,3] [,4] [,5]
##  [1,]    2    1    4    5    3
##  [2,]    2    1    5    3    4
##  [3,]    2    3    1    5    4
##  [4,]    2    3    4    5    1
##  [5,]    2    3    5    1    4
##  [6,]    2    4    1    5    3
##  [7,]    2    4    5    1    3
##  [8,]    2    4    5    3    1
##  [9,]    2    5    1    3    4
## [10,]    2    5    4    1    3
## [11,]    2    5    4    3    1
## [12,]    3    1    2    5    4

Generating Functions

A generating function is a power series in one indeterminate, whose coefficients encode information about a sequence of numbers \(a_n\) that is indexed by the natural numbers. They are an important tool in combinatorics.

Function \(G(x) = a_0 + a_1x + a_2x^2 + \ldots\) is the generating function of sequence \(\{a_0,a_1,a_2,\ldots\}\).

They can be used to awswer questions like this:

There are 3 pieces of 1 gram weight, 4 pieces of 2 grams weight, 2 pieces of 4 grams weight, how many different ways to weigh the value of 6-gram?

library(rSymPy) 

x <- Var("x")
sympy("a = 1+x+x**2+x**3")          # up to 3 pieces of 1 gram
## [1] "1 + x + x**2 + x**3"
sympy("b = 1+x**2+x**4+x**6+x**8")  # up to 4 pieces of 2 gram
## [1] "1 + x**2 + x**4 + x**6 + x**8"
sympy("c = 1+x**4+x**8")            # up to 2 pieces of 4 gram
## [1] "1 + x**4 + x**8"

sympy("simplify(a*b*c)")            # the generating function for this problem
## [1] "1 + x + 2*x**2 + 2*x**3 + 3*x**4 + 3*x**5 + 4*x**6 + 4*x**7 + 5*x**8 + 5*x**9 + 5*x**10 + 5*x**11 + 4*x**12 + 4*x**13 + 3*x**14 + 3*x**15 + 2*x**16 + 2*x**17 + x**18 + x**19"

The answer for the number of ways is in the coefficient of \(x^6\), in this case \(4\). The solutions are \(4+1+1; 2+4; 2+2+1+1; 2+2+2\).

We can extract these coefficients using R:

sol <- sympy("simplify(a*b*c)")            # the generating function for this problem
sol
## [1] "1 + x + 2*x**2 + 2*x**3 + 3*x**4 + 3*x**5 + 4*x**6 + 4*x**7 + 5*x**8 + 5*x**9 + 5*x**10 + 5*x**11 + 4*x**12 + 4*x**13 + 3*x**14 + 3*x**15 + 2*x**16 + 2*x**17 + x**18 + x**19"
monomials <- strsplit(sol,"+", fixed=TRUE)[[1]]
monomials[7]
## [1] " 4*x**6 "
regmatches(monomials[7],gregexpr("[0-9]+", monomials[7]))[[1]][1]
## [1] "4"

Another eg is to get the number of possible results from playing with several dice:

sympy("dice = x+x**2+x**3+x**4+x**5+x**6")
## [1] "x + x**2 + x**3 + x**4 + x**5 + x**6"
sympy("simplify(dice*dice)")   # two dice
## [1] "x**2 + 2*x**3 + 3*x**4 + 4*x**5 + 5*x**6 + 6*x**7 + 5*x**8 + 4*x**9 + 3*x**10 + 2*x**11 + x**12"
sympy("simplify(dice**3)")     # three dice
## [1] "x**3 + 3*x**4 + 6*x**5 + 10*x**6 + 15*x**7 + 21*x**8 + 25*x**9 + 27*x**10 + 27*x**11 + 25*x**12 + 21*x**13 + 15*x**14 + 10*x**15 + 6*x**16 + 3*x**17 + x**18"

or a coin eg:

We have three type of coins, 10 cents, 20 cents and 50 cents. What are the possible combinations for amounts up to 100 cents?

sympy("coin10 = 1+x**10+x**20+x**30+x**40+x**50+x**60+x**70+x**80+x**90+x**100")
## [1] "1 + x**10 + x**20 + x**30 + x**40 + x**50 + x**60 + x**70 + x**80 + x**90 + x**100"
sympy("coin20 = 1+x**20+x**40+x**60+x**80+x**100")
## [1] "1 + x**20 + x**40 + x**60 + x**80 + x**100"
sympy("coin50 = 1+x**50+x**100")
## [1] "1 + x**50 + x**100"
sympy("simplify(coin10*coin20*coin50)") # it gives more results, up to 300 cents
## [1] "1 + x**10 + 2*x**20 + 2*x**30 + 3*x**40 + 4*x**50 + 5*x**60 + 6*x**70 + 7*x**80 + 8*x**90 + 10*x**100 + 10*x**110 + 11*x**120 + 11*x**130 + 12*x**140 + 12*x**150 + 12*x**160 + 11*x**170 + 11*x**180 + 10*x**190 + 10*x**200 + 8*x**210 + 7*x**220 + 6*x**230 + 5*x**240 + 4*x**250 + 3*x**260 + 2*x**270 + 2*x**280 + x**290 + x**300"

Some useful relations:

Say I want to split \(n\) into summations of \(1,2,3,\ldots,m\). The generating function would be:

\(G(x) = (1+x+x^2+\ldots)(1+x^2+x^4+\ldots)\ldots(1+x^m+x^{2m}+\ldots) = (1-x)^{-1}(1-x^2)^{-1}\ldots(1-x^m)^{-1}\)

The generating function \(\frac{1}{(1-x)^2}\) generates what sequence?

sympy("ones = 1+x+x**2+x**3+x**4+x**5+x**6+x**7+x**8") # test first elements
## [1] "1 + x + x**2 + x**3 + x**4 + x**5 + x**6 + x**7 + x**8"
sympy("simplify(ones**2)") 
## [1] "1 + 2*x + 3*x**2 + 4*x**3 + 5*x**4 + 6*x**5 + 7*x**6 + 8*x**7 + 9*x**8 + 8*x**9 + 7*x**10 + 6*x**11 + 5*x**12 + 4*x**13 + 3*x**14 + 2*x**15 + x**16"

We can see that the sequence is \(\mathbb{N}\), ie, \(\{1,2,3,4,\ldots\}\)

Find the sequence \(\{a_n\}\) generated by \(\frac{3+78x}{1-3x-54x^2}\)

\[G(x) = \frac{3+78x}{1-3x-54x^2} = \frac{3+78x}{(1+6x)(1-9x)} = \frac{A}{1+6x}+\frac{B}{1-9x}\]

multiplying both fractions, then in the numerator we get a 2-degree polynomial that can be used to compare with \(3+78x\) to find both A and B:

\[G(x) = \frac{7}{1+6x}-\frac{4}{1-9x} = 7[1+(9x)+(9x)^2+\ldots]-4[1+(-6x)+(-6x)^2+\ldots]\]

This results in \(\{a_n\} = 7*9^n - 4*(-6)^n\)

This is the generation function of \(C(n,k)\).

Relation to Recurrence Relations

Ginve the recurrence relation \(h(n) = 2h(n-1) + 1\) with \(h(0)=0, h(1)=1\), how can we find the corresponding generating function given \(h(n)\)?

We want \(G(x) = h(1)x + h(2)x^2 + h(3)x^3 + \ldots\) (notice that \(h(0)=0\))

If we sum that with \(-2xG(x) = -2(h1)x^2 -2h(2)x^3 + \ldots\)

We get \((1-2x)G(x) = h(1)x - [h(2)-2h(1)]x^2 + [h(3)-2h(2)]x^3 + \ldots\)

Since \(h(n)-2h(n-1)=1\), \((1-2x)G(x) = x + x^2 + x^3 + \ldots = \frac{x}{1-x}\) and we have the expression for G(x).

Now, with \(G(x)\) we are able to find a closed expression for the sequence produced by \(h(x)\)

\[ \begin{array}{lclr} G(x) & = & \frac{x}{(1-x)(1-2x)} & \\ & = & \frac{1}{1-2x}-\frac{1}{1-x} & \\ & = & (1+2x+2^2x^2 + 2^3x^3 + \ldots) - (1+x+x^2+x^3+\ldots) & \\ & = & (2-1)x + (2^2-1)x^2 + (2^3-1)x^3 + \ldots & \\ & = & \sum_{k=1}^\infty (2^k-1)x^k & \\ \end{array} \]

So the sequence for \(h(n)\) is \(\{a_n\} = 2^n-1\)

Exponential Generating Function

For an exponential generating function \(G_e(x)\) we construct the following function:

\[G_e(x) = a_0 + \frac{a_1}{1!}x + \frac{a_2}{2!}x^2 + \frac{a_3}{3!}x^3 + \ldots\]

while generating functions are more used in combinatorial contexts, exponential generating functions are useful in contexts dealing with permutations.

An eg of use is to count permutations with repetitions, ie, multisets.

We want to know the possibilities of filling five number slots with number 1 to 4, considering that 1 must appear once or twice, 2 cannot appear more than once, 3 can appear up to three times, and 4 can only appear an even number of times.

We can compute the following \(G_e\):

sympy("number1 =     x/1 + x**2/2")
## [1] "x + x**2/2"
sympy("number2 = 1 + x/1")    
## [1] "1 + x"
sympy("number3 = 1 + x/1 + x**2/2 + x**3/6")
## [1] "1 + x + x**2/2 + x**3/6"
sympy("number4 = 1 +     + x**2/2 +         x**4/24")
## [1] "1 + x**2/2 + x**4/24"

sympy("simplify(number1*number2*number3*number4)") 
## [1] "x + 5*x**2/2 + 3*x**3 + 8*x**4/3 + 43*x**5/24 + 43*x**6/48 + 17*x**7/48 + 29*x**8/288 + x**9/48 + x**10/288"

The number we want is the coefficient of \(x^5\) which is \(43/24\). Since this is an exponential generating function, the \(x^5\) coefficient must be divided by \(5!=120\), so \(43x^5/24 = 215x^5/5!\). So, \(215\) is the number of ways to combine the numbers according to the constraints.

The exponential term is related to the Taylor expansion of \(e^x\):

\[e^x = \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \ldots\]

So, \(e^x\) generates \(\{1,1,1,\ldots\}\).

Some related expansions:

  • \(e^{-1}\) generates \(\{1,-1,1,-1,\ldots\}\)

  • \(e^{kx}\) generates \(\{1,k,k^2,k^3,\ldots\}\)

  • \(\frac{e^x+e^{-x}}{2}\) generates \(\{1,0,1,0,1,0,\ldots\}\)

Catalan Numbers

Catalan numbers appear in many counting problems. The n-th catalan number \(C_n, n\gt 0\) is

\[C_n = \frac{1}{n+1} {2n \choose n} = \frac{(2n)!}{(n+1)!n!} = \prod_{k=2}^n \frac{n+k}{k}\]

cn <- function(n) {
  choose(2*n,n) / (n+1)
}

xs <- 1:14
print(data.frame(n=xs, cn=cn(xs)), row.names=FALSE, quote=FALSE)
##   n      cn
##   1       1
##   2       2
##   3       5
##   4      14
##   5      42
##   6     132
##   7     429
##   8    1430
##   9    4862
##  10   16796
##  11   58786
##  12  208012
##  13  742900
##  14 2674440

Applications

A Dyck word is a string consisting of n X’s and n Y’s such that no initial segment of the string has more Y’s than X’s. An instance is the language of parenthesis with X=( and Y=).

There are six people in a library queuing up, three of them want to return the book “Interviewing Skills”, and 3 of them want to borrow the same book. If at the beginning, all the books of “Interviewing Skills” are out of stock in the library, how many ways can these people line up?

The result is based on \(C_3=5\) which we must multiply after considering the permutations of different people (three to deliver, three to borrow). So, \(C_3*3!*3! = 180\) ways to these people to line up.

To compute the solutions in R:

S <- c(1,2,3,-1,-2,-3)  # dyck words are represented as 'X' > 0 and 'Y' < 0
r <- length(S)

result <- permutations(r, r, S)  # get all possible r! permutations 
head(result)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]   -3   -2   -1    1    2    3
## [2,]   -3   -2   -1    1    3    2
## [3,]   -3   -2   -1    2    1    3
## [4,]   -3   -2   -1    2    3    1
## [5,]   -3   -2   -1    3    1    2
## [6,]   -3   -2   -1    3    2    1
tail(result)
##        [,1] [,2] [,3] [,4] [,5] [,6]
## [715,]    3    2    1   -3   -2   -1
## [716,]    3    2    1   -3   -1   -2
## [717,]    3    2    1   -2   -3   -1
## [718,]    3    2    1   -2   -1   -3
## [719,]    3    2    1   -1   -3   -2
## [720,]    3    2    1   -1   -2   -3

# filter only admissible results, ie, rows representing dyck words
#  1. convert negatives to -1, positives to +1
#  2. for each row, make a cumulative sum & find the min
#  3. if min is negative, the row represents an unbalanced expression & is filtered out
result <- result[apply(result, 1, function(x) min(cumsum(ifelse(x>0,1,-1))) ) >= 0,]
head(result)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1   -3    2   -2    3   -1
## [2,]    1   -3    2   -1    3   -2
## [3,]    1   -3    2    3   -2   -1
## [4,]    1   -3    2    3   -1   -2
## [5,]    1   -3    3   -2    2   -1
## [6,]    1   -3    3   -1    2   -2
tail(result)
##        [,1] [,2] [,3] [,4] [,5] [,6]
## [175,]    3    2    1   -3   -2   -1
## [176,]    3    2    1   -3   -1   -2
## [177,]    3    2    1   -2   -3   -1
## [178,]    3    2    1   -2   -1   -3
## [179,]    3    2    1   -1   -3   -2
## [180,]    3    2    1   -1   -2   -3

If we just want the Catalan sequences (X objects are identical, Y objects are identical) we can use the code for multisets and filter just the dyck words:

S <- c(1,1,1,-1,-1,-1)  # dyck words are represented as 'X'=+1 and 'Y'=-1
r <- length(S)

result <- unique(permutations(r, r, S, set=FALSE))  # get just the different permutations 
head(result,10)
##       [,1] [,2] [,3] [,4] [,5] [,6]
##  [1,]    1    1    1   -1   -1   -1
##  [2,]    1    1   -1    1   -1   -1
##  [3,]    1    1   -1   -1    1   -1
##  [4,]    1    1   -1   -1   -1    1
##  [5,]    1   -1    1    1   -1   -1
##  [6,]    1   -1    1   -1    1   -1
##  [7,]    1   -1    1   -1   -1    1
##  [8,]    1   -1   -1    1    1   -1
##  [9,]    1   -1   -1    1   -1    1
## [10,]    1   -1   -1   -1    1    1
result <- result[apply(result, 1, function(x) min(cumsum(x)) ) >= 0,]
head(result)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    1   -1   -1   -1
## [2,]    1    1   -1    1   -1   -1
## [3,]    1    1   -1   -1    1   -1
## [4,]    1   -1    1    1   -1   -1
## [5,]    1   -1    1   -1    1   -1

Stirling Numbers

First Kind

Sitrling numbers also appear in many application of combinatorics.

The signed Stirling numbers of the first kind \(s(n,m)\) are defined such that the number of permutations of n elements which contain exactly \(m\) permutation cycles is the nonnegative number \(|s(n,m)|\) where,

\[s(n+1,m)=s(n,m-1)-n s(n,m)\]

with \(m \gt 0\), and \(s(0,0)=1\), \(s(n,0)=s(0,n)=0, n \gt 0\)

In R:

library(copula)

stirs1 <- matrix(rep(NA,100), ncol=10)
for(n in 0:9) {
  for(m in 0:n) 
     stirs1[n+1,m+1] <- abs(Stirling1(n,m))
}
print(stirs1, na.print=" " )
##       [,1]  [,2]   [,3]   [,4]  [,5]  [,6] [,7] [,8] [,9] [,10]
##  [1,]    1                                                     
##  [2,]    0     1                                               
##  [3,]    0     1      1                                        
##  [4,]    0     2      3      1                                 
##  [5,]    0     6     11      6     1                           
##  [6,]    0    24     50     35    10     1                     
##  [7,]    0   120    274    225    85    15    1                
##  [8,]    0   720   1764   1624   735   175   21    1           
##  [9,]    0  5040  13068  13132  6769  1960  322   28    1      
## [10,]    0 40320 109584 118124 67284 22449 4536  546   36     1

notice that R arrays start at index 1. So \(s(7,1)=24\) is at position \([6][2]\).

There is a formula to compute \(s(n,2)\) that uses harmonic numbers. The n-th harmonic number, \(H_n\), is the sum of the reciprocals of the first n natural numbers, ie,

\[H_n = \sum_{i=1}{n} \frac{1}{n}\]

the formula is

\[s(n,2) = n! H_n\]

Second Kind

A Stirling number of the second kind (or Stirling partition number), \(S(n,m)\) is the number of ways to partition a set of \(n\) objects into \(m\) non-empty subsets.

\[S(n,m) = \frac{1}{m!} \sum_{j=0}^m (-1)^j {m \choose j} (m-j)^n\]

or following the recurrence relation

\[S(n.m) = S(n-1,m-1) + n S(n-1,m)\]

with \(S(0,0)=1\) and \(S(n,0)=0, n \gt 0\)

In R:

stirs2 <- matrix(rep(NA,100), ncol=10)
for(n in 0:9) {
  for(m in 0:n) 
     stirs2[n+1,m+1] <- Stirling2(n,m)
}
print(stirs2, na.print=" " )
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1                                              
##  [2,]    0    1                                         
##  [3,]    0    1    1                                    
##  [4,]    0    1    3    1                               
##  [5,]    0    1    7    6    1                          
##  [6,]    0    1   15   25   10    1                     
##  [7,]    0    1   31   90   65   15    1                
##  [8,]    0    1   63  301  350  140   21    1           
##  [9,]    0    1  127  966 1701 1050  266   28    1      
## [10,]    0    1  255 3025 7770 6951 2646  462   36     1

How many ways can you divide four people A B C and D into two groups?

Stirling2(4,2)
## [1] 7

The solutions are {A,B},{C,D};{A,C},{B,D};{A,D},{B,C};{A},{B,C,D};{B},{A,C,D};{C},{A,B,D};{D},{A,B,C}

Other formulas:

  • The ways to place n distict balls into m distinct boxes (leaving no empty box) is \(m! S(n,m)\).

  • The ways to place n distict balls into m identical boxes (leaving no empty box) is \(S(n,m)\).

  • The ways to place n distict balls into m distinct boxes (empty boxes allowed) is \(m^n\).

  • The ways to place n distict balls into m identical boxes (empty boxes allowed) is \(S(n,1)+S(n,2)+\ldots+S(n,k)\) where \(k=m \gt n ? n : m\).

  • The ways to place n identical balls into m distinct boxes (leaving no empty box) is \(C(n-1,m-1)\).

  • The ways to place n identical balls into m identical boxes (leaving no empty box) is is the coefficient of \(x^n\) in \(G(x) = \frac{x^m}{(1-x)(1-x^2)\cdots(1-x^m)}\)..

  • The ways to place n identical balls into m distinct boxes (empty boxes allowed) is \(C(n+m-1,n)\).

  • The ways to place n identical balls into m identical boxes (empty boxes allowed) is the coefficient of \(x^n\) in \(G(x) = \frac{1}{(1-x)(1-x^2)\cdots(1-x^m)}\).

Bell Numbers

Bell numbers follow the recurrence relation

\[B_n = \sum_{i=1}^n {n-1 \choose i-1} B_{n-i}\]

with \(B_0 = 1\) (only one partition of the empty set).

# returns a vector with the first (max+1) bell numbers
# pre: max > 0
bell <- function(max) {
  result <- c(1,rep(NA,max-1)) # B(0) = 1
  
  for (n in 1:max) {
    is = 1:n
    result[n+1] <- sum( choose(n-1,is-1)*result[1+n-is] )
  }
  
  result
}

print(data.frame(n=0:12, Bn=bell(12)), digits=7, row.names=FALSE)
##   n      Bn
##   0       1
##   1       1
##   2       2
##   3       5
##   4      15
##   5      52
##   6     203
##   7     877
##   8    4140
##   9   21147
##  10  115975
##  11  678570
##  12 4213597

The Bell number \(B_n\) gives the number of partitions of a set of size \(n\).

Eg, for set \(\{1,2,3\}\), with \(n=3\), \(B_3=5\) since there are five different ways to partition it (the order is irrelevant),

\[\{ \{a\}, \{b\}, \{c\} \}~;~ \{ \{a\}, \{b, c\} \}~;~ \{ \{b\}, \{a, c\} \}~;~ \{ \{c\}, \{a, b\} \}~;~ \{ \{a, b, c\} \}\]

bell(3)[4] # returns B_3
## [1] 5

They can also be computed via Stirling numbers of the 2nd kind:

\[B_n = \sum_{k=0}^n S(n,k)\]

library(partitions) # http://cran.r-project.org/web/packages/partitions/vignettes/setpartitions.pdf

listParts(3)
## [[1]]
## [1] (1,2,3)
## 
## [[2]]
## [1] (1,3)(2)
## 
## [[3]]
## [1] (1,2)(3)
## 
## [[4]]
## [1] (2,3)(1)
## 
## [[5]]
## [1] (1)(2)(3)
# returns, for each row, from which partition a given element belongs to
as.matrix(t(setparts(3)), ncol=3)
##      [,1] [,2] [,3]
## [1,]    1    1    1
## [2,]    1    2    1
## [3,]    1    1    2
## [4,]    2    1    1
## [5,]    1    2    3

Partitions of an Integer

The five partitions of \(4\) are \(4=3+1=2+2=2+1+1=1+1+1+1\).

The number of partitions of an integer is given by the partition function p(n) which does not have a closed form.

However we can compute its value using the following recurrence relation

\[p(n) = \sum_{k \in \mathbb{Z}-\{0\}} (-1)^{k-1} p\big(n-\frac{k(3k-1)}{2}\big) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(n-12) \ldots \]

with \(p(0)=1\) and \(p(n)=0, n \lt 0\).

Package partitions can used to enumerate the partitions of an integer:

library(partitions)

n <- 4
m <- as.matrix(t(parts(n)), ncol=n)
m[m==0] = NA
print(m, na.print="")
##      [,1] [,2] [,3] [,4]
## [1,]    4               
## [2,]    3    1          
## [3,]    2    2          
## [4,]    2    1    1     
## [5,]    1    1    1    1
P(4) # just the number of partitions of given integer
## [1] 5

If the ordering matters (ie, \(1+4 \neq 4+1\)):

t(compositions(5))
##                
##  [1,] 5 0 0 0 0
##  [2,] 1 4 0 0 0
##  [3,] 2 3 0 0 0
##  [4,] 1 1 3 0 0
##  [5,] 3 2 0 0 0
##  [6,] 1 2 2 0 0
##  [7,] 2 1 2 0 0
##  [8,] 1 1 1 2 0
##  [9,] 4 1 0 0 0
## [10,] 1 3 1 0 0
## [11,] 2 2 1 0 0
## [12,] 1 1 2 1 0
## [13,] 3 1 1 0 0
## [14,] 1 2 1 1 0
## [15,] 2 1 1 1 0
## [16,] 1 1 1 1 1

The package is also able to return the partitions restricted to a certain maximum size:

t(restrictedparts(6,3))    # partition of 6 up to 3 partitions
##           
## [1,] 6 0 0
## [2,] 5 1 0
## [3,] 4 2 0
## [4,] 3 3 0
## [5,] 4 1 1
## [6,] 3 2 1
## [7,] 2 2 2
R(3,6, include.zero =TRUE) # just the number
## [1] 7

Permutation Group

A permutation group on set \(1, 2, \ldots, n\) is a 1-1 mapping on itself. It is represented by \(\left( \begin{smallmatrix} 1 & 2 & \ldots & n \cr a_1 & a_2 & \ldots & a_n \end{smallmatrix} \right)\) where \(a_1a_2\ldots a_n\) is a set arragement. The same permutation may have \(n!\) representations. Eg, these two representations are equivalent:

\[\left( \begin{matrix} 1 & 2 & 3 & 4 \cr 3 & 1 & 2 & 4 \end{matrix} \right) = \left( \begin{matrix} 3 & 1 & 4 & 2 \cr 2 & 3 & 4 & 1 \end{matrix} \right) \]

Permutation multiplication receives two permutation groups and return the two mappings composition.

Eg:

\[\left( \begin{matrix} 1 & 2 & 3 & 4 \cr 3 & 1 & 2 & 4 \end{matrix} \right) . \left( \begin{matrix} 3 & 1 & 2 & 4 \cr 2 & 4 & 3 & 1 \end{matrix} \right) = \left( \begin{matrix} 1 & 2 & 3 & 4 \cr 2 & 4 & 3 & 1 \end{matrix} \right) \]

This operation is not commutative but it defines a group with the permutations of a given set.

A quick reminder: given a set \(G\) and an operation \((.) : G \times G \rightarrow G\), the tuple \((G,.)\) is a group if it has the next four properties:

The identity element is \[e= \left( \begin{matrix} 1 & 2 & 3 & 4 \cr 1 & 2 & 3 & 4 \end{matrix} \right)\] and the inverse is \[\left( \begin{matrix} 1 & 2 & 3 & 4 \cr a_1 & a_2 & \ldots & a_n \end{matrix} \right)^{-1} = \left( \begin{matrix} a_1 & a_2 & \ldots & a_n \cr 1 & 2 & 3 & 4 \end{matrix} \right) \]

The group of all permutations of \(\{1,2,\ldots,n\}\) over the composition is known as n-order symmetric group, or just \(S_n\).

A more compact notation for, say, \(e= \left( \begin{matrix} 1 & 2 & 3 & 4 & 5\cr 5 & 2 & 3 & 1 & 4 \end{matrix} \right)\) is \((154)(2)(3)\) which is made by starting at the first position, \(1\), and travel the sequence \(1 \rightarrow 5 \rightarrow 4 \rightarrow 1\) until we reach a cycle; then repeat for the elements not included yet.

A theorem states that given \(p=(a_1a_2\cdots a_n)\), \(p^n=(1)(2)\cdots(n)=e\).

Eg: \(p=(123)\), then \(p^2=p.p= \left( \begin{matrix} 1 & 2 & 3\cr 2 & 3 & 1 \end{matrix} \right) . \left( \begin{matrix} 1 & 2 & 3\cr 2 & 3 & 1 \end{matrix} \right) = \left( \begin{matrix} 1 & 2 & 3\cr 3 & 1 & 2 \end{matrix} \right) =(132)\). And \(p^3=p.p^2= \left( \begin{matrix} 1 & 2 & 3\cr 1 & 2 & 3 \end{matrix} \right) = (1)(2)(3)=e\)

Def: for a \(p \in S_m\), the cyclic group generated by p is \(\left<p\right> = \{ p^n | 0 \leq n \lt k, p^k=e \}\)

Theorem: For \(p\in S_m\), \(\left<p\right>\) is a subgroup of \(S_m\).

Conjugacy Class

Every permutation has format

\[(a_1a_2\cdots a_{k_1}) (b_1b_2\cdots b_{k_2}) \cdots (h_1h_2\cdots h_{k_m})\]

where \(\sum_i k_i = n\).

For \(S_n\) we can use a more compact notation,

\[1^{c_1} 2^{c_2} \cdots n^{c_n}\]

where \(\sum_k k.c_k = n\)

Eg, in \(S_5\), \(1+1+3\) can be represented by \(1^23^1\). The five partitions of \(S_5\) would be written as \(4^1; 3^11^1; 2^2, 2^11^2, 1^4\).

The conjugacy class over group \(G\) of element \(a\in G\) is

\[Cl(a) = \{ g\in G | x\in G, g = xax^{-1} \}\]

All pair of elements \(a,b\in G\) either belong to the same equivalence class and then \(Cl(a)=Cl(b)\), or do not and then \(Cl(a) \cap Cl(b)=\emptyset\).

A conjugate class is an equivalent relation, partitioning \(G\) into equivalence classes.

Another reminder: A relation \(R\) on set \(S\) is an equivalence relation iff

An equivalence class of \(x\in S\) is defined as \(\{y \in S | yRx \}\).

For symmetric groups, permutations with the same notation \(1^{c_1} 2^{c_2} \cdots n^{c_n}\) belong to the same conjugacy class.

Egs:

Symmetric group \(S_3\) has three conjugacy classes: \(1^3; 1^12^1; 3^1\).

Symmetric group \(S_4\) has five conjugacy classes: \(1^4; 1^22^1; 2^2; 1^13^1; 4^1\).

The conjugacy class \(2^2\) of \(S_4\) is composed of 3 different permutations: \((12)(34); (13)(24); (14)(23)\).

The conjugacy class \(1^13^1\) of \(S_4\) is composed of 8 different permutations: \((123)(4); (124)(3); (132)(4); (134)(2); (142)(3); (143)(2); (234)(1); (243)(1)\).

To compute the conjugacy class size \(1^{c_1} 2^{c_2} \cdots n^{c_n}\) of the symmetric group \(S_n\) use the following formula:

\[\frac{n!}{\prod_{k=1}^n k^{c_k} c_k!}\]

Eg: conjugacy class \(1^13^1\) of \(S_4\) has \(\frac{4!}{1^1\times 1! \times 3^1 \times 1!} = 24/3 = 8\) elements.

Eg2: conjugacy class \(1^22^24^1\) of \(S_{10}\) has \(\frac{10!}{1^2 \times 2! \times 2^2 \times 2! \times 4^1 \times 1!}=56700\) elements.

The number of permutation in \(S_n\) with \(m\) cycles is given by the Stirling number of the first kind, \(s(n,m)\).

Eg: \(s(5,4)=10\) which corresponds to the 4-cycle permutations on \(S_5\), namely \((12)(3)(4)(5); (2)(13)(4)(5); (1)(2)(3)(14)(5); \ldots (1)(2)(3)(45)\).

Stabilizers & Orbits

Let \(G\) be a permutation group on set \(S\), and \(x\in S\). The set

\[G_x = {g \in G | g(x) = x}\]

is called the stabilizer of \(x\).

Eg: If \(G = \{ e, (12)(3)(4), (1)(2)(34), (12)(34) \}\), \(G_1 = \{ e, (1)(2)(34)\}\)

The subset of all images of \(x\in S\) under permutations of group \(G\), is called the orbit of \(x\) in \(G\),

\[G(x) = \{ g(x) | g \in G \}\]

Orbits are equivalence classes of \(G\), and follow the property \(|G| = |G_x| |G(x)|\) known as the Orbit-stabilizer theorem.

Eg: If \(G = \{ e, (12)(3)(4), (1)(2)(34), (12)(34) \}\), \(G(1) = \{ 1, 2\}\) since for \(g=(12)(3)(4), g(1)=2\) and for \(g=e, g(1)=1\) (no other values are possible for this \(G\)).

The previous property also holds: \(|G|=4\), \(|G_1|\times|G(1)|=4\)

Inclusion-Exclusion Principle

There are two fundamental counting principles

However if \(A\) and \(B\) are not mutually exclusive we cannot use these principles.

The Inclusion-Exclusion Principle states that if \(A \cap B \neq \emptyset\) where the intersection has size \(k\), then event \(A \lor B\) can happen \(m+n-k\) ways, ie, \(|A \cup B| = |A| + |B| - |A \cap B|\).

For three sets we’ll have \(|A \cup B \cup C| = |A|+|B|+|C|-|A \cap B|-|A \cap C|-|A \cap B|+|A\cap B \cap C|\).

Calculate the size of S, the permutation set of which do not contain “ace” and “df”.

The permutations of six letters is \(6!=720\)

Assume A is the permutation set which “ace” is an element. If we consider it as one element, the permutation set only has four elements, so \(|A|=4!\). Assume B is the permutation set which “df” is an element, so \(|B|=5!\). Also consider \(|A \cap B|\) which is the set with both “ace” and “df”. Using the same reasoning, \(|A \cap B|=3!\).

So, \(|S| = 6! - |A| - |B| + |A \cap B| = 720 - 24 - 120 + 6= 582\) elements.

How many primes are smaller than 120?

Since \(11^2 = 121\) the composite numbers always have a divisor smaller than \(11\).

Say \(A_i\) is the set of number divisible by \(i\) smaller than 120, where \(i=2,3,5,7\).

\[|A_i| = \lfloor\frac{120}{i}\rfloor\], so \(|A_2|=60, |A_3|=40, |A_5|=24, |A_7|=17\).

We must check the intersections, say \(A_2\cap A_3\) is the set of number divisible by \(6\), so

\[|A_2\cap A_3| = \lfloor\frac{120}{2\times 3}\rfloor = 20\]

and \(|A_2\cap A_5| = 12\), \(|A_2\cap A_7| = 8\), \(|A_3\cap A_5| = 8\), \(|A_3\cap A_7| = 5\), \(|A_5\cap A_7| = 3\).

Now, we need to (re)add the sets \(|A_2\cap A_3\cap A_5|=4\), \(|A_2\cap A_3\cap A_7|=2\), \(|A_2\cap A_5\cap A_7|=1\) and \(|A_3\cap A_5\cap A_7|=1\)

And discount again \(|A_2\cap A_3\cap A_5\cap A_7|=0\).

After all the summations and subtractions, we get \(27\) primes.

A useful function related to this problem is the Euler function \(\Phi(n)\) is calculates the number of integers smaller than \(n\) and relatively prime to \(n\):

\[\Phi(n) = n(1-\frac{1}{p_1})(1-\frac{1}{p_2})\cdots(1-\frac{1}{p_k1})\]

where \(p_i\) is the i-th prime (\(p_1=2, p_2=3, \ldots\)) that are included in the prime expansion of \(n\).

Eg: To calculate \(\Phi(60)\), first decompose the number \(60 = 2^2 \times 3 \times 5\), then \(\Phi(60) = 60(1-1/2)(1-1/3)(1-1/5) = 16\)

Pigeonhole Principle

The Pigeonhole Principle states that if \(n\) items are each put into one of \(m\) boxes, where \(n>m\), then one of the boxes contains more than one item.

More generally, if we have \(m\) boxes and \(kn+1\) items, there must be a box with at least \(k+1\) items.

Ramsey Numbers

An eg of Ramsey problem: among 6 people there must be 3 mutual friends or 3 mutual strangers. This can be seen as a coloring problem of complete graphs:


Paint each edge red or blue. There will always be a monochromatic triangle.

Ramsey’s theorem states that one will find monochromatic cliques in any edge colouring of a sufficiently large complete graph. To demonstrate the theorem for two colours (say, blue and red), let \(r\) and \(s\) be any two positive integers. Ramsey’s theorem states that there exists a least positive integer \(R(r,s)\) for which every blue-red edge colouring of the complete graph on \(R(r,s)\) vertices contains a blue clique on \(r\) vertices or a red clique on \(s\) vertices. Here \(R(r,s)\) signifies an integer that depends on both \(r\) and \(s\). ref

The previous example is written \(R(3,3)=6\). These numbers are called Ramsey Numbers.

Another eg is \(R(4,2)=4\) meaning that in a complet four node graph (called \(K_4\)) there is a clique of size four or a clique of size two.

The first complete graphs:

This is a very complex problem, we only know a small set of Ramsey numbers. Even \(R(5,5)\) is unknown, but we know that it must lie between \([43,49]\). Detailed info at Wolfram’s Mathworld.