Implementing Special Numbers: Advanced List Manipulations

Keehong Song
Pusan National University,South Korea
khsong@web.pusan.ac.kr

This is a Mathematica Notebook which implements some special recurrence relations in combinatorics. It provides additional list manipulation capabilities to complement the Mathematica kernel and 'DiscreteMath' standard package. Also in the Notebook, the utilties of the routines are demontrated as examples. The first two functions, 'stirling2nd' and 'Bell' are about the partition of a set while the rest about the permutation.

References:
Graham,R., Knuth D., and Patashnik, O., Concrete Mathematics, Addison-Wesley,1990.
Grimaldi,R.Discrete and Combinatorial Mathematics,3rd ed.,Addison-Wesley,1994.

Definition (Stirling number of the second kind)
The number of partitions of the n objects into k nonempty subsets is denoted by{({{n}, {k}})}.
Then the recurrence relation holds
	{({{n}, {k}})}={({{n - 1}, {k - 1}})}+k{({{n - 1}, {k}})}, {({{n}, {1}})}={({{n}, {n}})}=1

stirling2nd[{}, k_Integer] := {{}} /; Length[l] == 0
stirling2nd[l_List, 1] := {{l}}
\(stirling2nd[l_List, 
        k_Integer] := 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {Partition[l, 1]} /; k == Length[l]\)
stirling2nd[l_List, k_Integer] := {{}} /; k > Length[l]
\(stirling2nd[l_List, k_Integer] := 
	
      Module[{i}, 
        Join[Flatten[
				\((Table[
                    Insert[#1, First[l], {i, 1}], {i, k}] &)\)
			 /@ \ 
              stirling2nd[Rest[l], k], 
            1], 
			\((Prepend[#1, {First[l]}] &)\) /@ \ 
            stirling2nd[Rest[l], k - 1]]]\)

Partition the set {a, b, c, d, e} into three nonempty subsets.

size3 = stirling2nd[{a, b, c, d, e}, 3]

Length[%]

The number of ways of partitioning can also be found using the built-in function, 'StirlingS2'.

StirlingS2[5, 3]

Example: Onto functions
For A = {a, b, c, d, e} and B = {1, 2, 3}, there are 150 onto functions from A to B. To obtain the number, many different ways are possible. Using the idea of partition, the set A can be partitioned into 3 nonempty sets and then these partitioned subsets are mapped one-to-one onto the codomain B.  So the number of onto function from A to B is  StirlingS2[5,3] × 3!=150

3! StirlingS2[5, 3]

Once subsets are indexed, then it is easy to see that the index itself represents each element in the codomain.

indexed = MapIndexed[List, #] & /@  size3;

We get all the functions by permuting the subscripts in the output.

Map[Subscripted[#1[Sequence @@ #2]] & @@ # &, 
      indexed, {2}]

Now we generate all the onto functions. To the end, first we define the routine, 'mapping', which gives actual ordered pairs.

\(mapping[domain_List, codomain_List]\  := 
	
      Flatten[
	\ MapThread[\ 
          Distribute[{#1, #2}, \ List] &, {domain, codomain}], 1]\)

The number of one-to-one functions from the indexed partition number {1, 2, 3} to elements in the codomain, {1, 2, 3}, is 3!

perm = Permutations[{1, 2, 3}]

Here we have onto functions we were after.

onto =  
Table[mapping[#, 
                Transpose[{perm[[i]]}]] &  /@  size3,  {i, 6 }]  // 
        Flatten[#, 1] &

Using the nonempty partition, we can develop another way of finding onto function. This time we use patterns and find the nonempty sets of the pattern,{{*},{*},{*,*}}.

Cases[size3, {{_}, {_}, {_, _, _}}]

The total number of the pattern is  ((5!)/(1! 1! 3!)).

Multinomial[1, 1, 3]

Now we find all the patterns of partitions possible.

The total number of the pattern is (( 
                5!)/(1! 2! 2!)).

Cases[size3, {{_}, {_, _}, {_, _}}]

The total number of the pattern is (( 
                5!)/(1! 3! 1!)).

Cases[size3, {{_}, {_, _, _}, {_}}]

The total number of the pattern is((5!)/(2! 1! 2!)).

Cases[size3, {{_, _}, {_}, {_, _}}]

The total number of the pattern is (( 
                5!)/(2! 2! 1!)).

Cases[size3, {{_, _}, {_, _}, {_}}]

The total number of the pattern is (( 
                5!)/(3! 1! 1!)).

Cases[size3, {{_, _, _}, {_}, {_}}]

Now we digress for a moment to get the solution to the equation, (Subscript[k, 1]) + (Subscript[k, 2]) + ⋯ + (Subscript[k, m]) = n, (Subscript[k, i]) ≥ 1, 
      i = 1  ⋯n.  Many elementary text books on the subject expains the way of finding the solution to the equation as follows:
for the case, n = 5, insert '+' sign in between the list of 1's, {1, 1, 1, 1, 1} and then perform addition of the group of 1's separated by the '+' signs.

ones = Table[1, {5}]

<< DiscreteMath`

The possible positions for '1' s to be inserted are {2, 3, 4, 5}.

pos = KSubsets[{2, 3, 4, 5}, 2]

Insert '1's into the set at the positions above.

s = Insert[ones, "+", Transpose[{#}]] &  /@  pos

Now perform addition of the 1's separated by the '+' marks.

Cases[
      s,  {a__, "+", b__, "+", c__}  :>  {Plus[a],  Plus[b],  
          Plus[c]}]

So far, we have demonstrated the concept behind the algorithm of finding the solution to the equation.

Partition[Join[{0}, #, {5}], 2, 1] &  /@  pos

Now we compute the difference of the elements in each pair in constructing a single function to find the solution set to the equation we started with.

\(toInteger[n_Integer, k_Integer]\  := 
	
      Map[
		\((Map[\((\(#2 - #1 &\) @@ #)\) &, \ 
			Partition[
                Join[{0}, #, {n}], 2, 1]\ ])\) &, 
		KSubsets[
          Range[n - 1], k - 1]
	]\)

toInteger[5, 3]

\(Now\ we\ are\ in\ a\ position\ to\ evaluate\ the\ number\ of\ onto\ functions\), (∑_({{(Subscript[k, i]) ≥ 1, i = 1  ⋯n}, {(Subscript[k, 1]) + (Subscript[k, 2]) + ⋯ + (Subscript[k, m]) = n}}))\(\(\(n!\)\/\(\(k\_1!\) \(k\_2!\) \(⋯k\_m!\)\)\)\(.\)\)

Multinomial @@ # &  /@  %

Plus @@ %

To verify validity of the method of counting.

3! StirlingS2[5, 3]

Definition (Bell number)
\(The\ number\ of\ partitions\ of\ a\ set\ with\ \ n\ element\ is\ denoted\ by\ \(B\_\(\(n\)\(\ \ \)\)\) 
        which\ satisfy\ the\ following\ recurrence\ relations\)
	 (Subscript[B, n])=(∑_(i = 0)^(n - 1))(({{n - 1}, {i}}))(Subscript[B, i])	
and in closed form, (Subscript[B, n]) is \(\(1\/e\) \(∑\+\(i = 0\)\%∞ i\^n\/\(i!\)\)\)

,which is in turn the same as (∑_({{(Subscript[k, i]) ≥ 0, i = 1  ⋯n}, {(Subscript[k, 1]) + 2  (Subscript[k, 2]) + ⋯ + (Subscript[nk, n]) = n}}))\(\(n!\)\/\(\(k\_1!\) \(k\_2!\) \(⋯k\_n!\) \(\((\(1!\))\)\^k\_1\) \(\((\(2!\))\)\^k\_2\) \(⋯(\(n!\))\)\^k\_n\)\)

<< DiscreteMath`
BellPartition[{}] := {{}}
BellPartition[{a_}] := {{{a}}}
\(BellPartition[l_List] := 
      Module[{i, 
          n = Length[l]}, 
		\((FlattenAt[#1, 2] &)\) /@ 
			Flatten[
            Table[\((Distribute[{{Flatten[{First[
                                l], #1}]}, 
			BellPartition[
                          Complement[Rest[l], #1]]}, 
                      List] &)\) /@ 
				KSubsets[Rest[l], i], {i, 0, 
                n - 1}], 2]]\)

The Bell number of the set {a, b, c, d} is the same as the number of the equivalence classes of the set.

class = BellPartition[{1, 2, 3, 4}]

Length[%]

Use closed form to find the number of equivalence classes.

bellCount[
        n_] := (1/E) (∑_(i = 0)^∞) ((i^n)/(i!))

bellCount[4]

And  we can get the number using the formula,(∑_({{(Subscript[k, i]) ≥ 0, i = 1  ⋯n}, {(Subscript[k, 1]) + 2  (Subscript[k, 2]) + ⋯ + n (Subscript[k, n]) = n}}))\(\(n!\)\/\(\(k\_1!\) \(k\_2!\) \(⋯k\_n!\) \(\((\(1!\))\)\^k\_1\) \(\((\(2!\))\)\^k\_2\) \(⋯(\(n!\))\)\^k\_n\)\), 
which is more intimidating  look .   The not so trivial part is to evaluate the formula .  'Discretemath' package function'Partitions' finds the solution set of the equation(Subscript[k, 1]) + 2  (Subscript[k, 2]) + ⋯ + n (Subscript[k, n]) = n. 'Partitions[n, k]' gives the partitions of the integer,\(n\ with\ the\ numbers\ less\ than\ or\ equal\ to\ \(\(k\)\(.\)\)\)

Partitions[4, 4]

freq = Table[Count[#, i], {i, 4}] & /@ %

\(\(4!\)/\(Product[\(#[\([i]\)]!\)\ \(i!\)\^#[\([i]\)], {i, 
              Length[#]}] &\) /@ freq\)

Plus @@ %

Once we get the partition of a set, obtaining the relations of the equivalence classes is simply a matter of performing the cartesian product of partitions.

pairs  =  
      Table[mapping[class[[i]], class[[i]]], {i, 3}]

The claim is easily checked as follows:

relations  =  FromOrderedPairs[#] &  /@  pairs

EquivalenceRelationQ  /@  relations

Now we retrieve the classes from the set of pairs.

Names["DiscreteMath`Combinatorica`*Q"]

EquivalenceClasses  /@  relations

Definition (Stirling number  of the first kind)
The number of permutations of the n objects having exactly k cycles is denoted by  [({{n}, {k}})].
	[({{n}, {k}})]=[({{n - 1}, {k - 1}})]+(k - 1)[({{n - 1}, {k}})], [({{n}, {0}})]=0,[({{n}, {n}})]=1

Stirling1st[{}, k_Integer] := {{}}
Stirling1st[l_List, 0] := {}
Stirling1st[l_List, k_Integer] := {Partition[l, 1]} /; 
        k == Length[l]
Stirling1st[l_List, k_Integer] := {{}} /; k > Length[l]
\(Stirling1st[l_List, k_Integer] := 
	
      Module[{i}, 
		Join[
          Flatten[\((
			Table[
                    Insert[#1, First[l], {i, j}], {i, k}, 
			{j, 
                      Length[#1[[
                          i]]]}] &)\) /@ 
              Stirling1st[Rest[l], k], 
            2], 
			\((Prepend[#1, {First[l]}] &)\) /@ 
            Stirling1st[Rest[l], k - 1]]]\)

As an example, the set of cycles of size 2 from the permutations of {1, 2, 3, 4}

cyc2  =  Stirling1st[{1, 2, 3, 4}, 2]

Length[%]

StirlingS1[4, 2]

Of course, the cycles of size 2 can be found from the total space of permutation of {1, 2, 3, 4}.

<< DiscreteMath`

allcyc = ToCycles /@ Permutations[{1, 2, 3, 4}]

Length[Select[allcyc, Length[#] == 2 &]]

Definition (derangements)
\(A\ derangement, 
        denoted\ by\ D\_n, \ \(is\ a\ permutation\ π\ with\ no\ fixed\ point;\)\) 
\(that\ is\ \(π \((x)\)\) = \(x\ has\ no\ solution . \ \ 

              Then\ D\_n = \((n - 1)\) \((D\_\(n - 1\) + 
                D\_\(n - 2\))\)\ holds . \ 
              The\ derangement\ is\ also\ called\ ' 
            problem\ de\ Recontre . '\)\)

recontre[{a_}] := {}
recontre[{a_, b_}] := {{b, a}}
\(recontre[l_List] := 
      Module[{i, n = Length[l]}, 
		Join @@ 
          Table[Join[
			Map[
                Prepend[#, 
                    l[[
                      i]]] &, 
			Map[
                  Insert[#, First[l], {i - 1}] &, 
				\ recontre[
                    Delete[l, {{1}, {i}}]]]], 
					Map[
                Prepend[#, 
                    l[[
                      i]]] &, 
				\ recontre[
                  Insert[Delete[l, {{1}, {i}}], 
					First[
                      l], {i - 1}]]]], 
								{i, 2, 
              n}]
				]\)

displaced = recontre[{1, 2, 3, 4}]

From all the permutations of the set {1, 2, 3, 4}, select only the ones with no fixed point.' DerangementQ' is a 'DiscreteMath'package function.

Select[Permutations[{1, 2, 3, 4}],  DerangementQ]

MapThread[#1 =!= #2 &, {#, Range[Length[#]]}] &  /@  
      Permutations[{1, 2, 3, 4}]

And @@ MapThread[#1 =!= #2 &, {#, Range[Length[#]]}] &  /@  
      displaced

We are certain that the results by the routine 'recontre' are correct.

And @@ %

We can also use generating function to find the derangements.  Define an indicator function X = (1 - (Subscript[X, 1])) (2 - 
              (Subscript[X, 2])) ⋯(n - (Subscript[X, n])). Then the function is nonzero only if each of  (Subscript[X, 1]), (Subscript[X, 2]), ⋯, (Subscript[X, n]) assumes different value from its own index.

indicator = Product[i - x[i], {i, 1, 5}]

perm = Permutations[{1, 2, 3, 4, 5}];

(Subscript[X, 1]), (Subscript[X, 2]), ⋯, (Subscript[X, n]) assumes the each element of the permutations.

vars = Array[x, {5}]

tORf = indicator /. Thread[vars -> #] & /@ perm

Length[Select[tORf, # ≠ 0 &]]

When just the number of derangements is needed, there is a simple formula for that.

(Subscript[D, n])= {({{⌊((n!)/e) + 
                          (1/2)⌋     , 
                      n ≳ 1}, {1                          , 
                      n = 0}})

numberOfRecontre[0] := 1
numberOfRecontre[k_] := 
      Floor[N[((k!)/E)] + (1/2)] /; k ≥ 1

numberOfRecontre[5]

Define h(n, k) to be the number of permutations, each of  which has k fixed points.
Then
	h(n, k)=(({{n}, {k}}))(Subscript[D, n - k])
holds

h[n_, k_] := Binomial[n, k] numberOfRecontre[n - k]

Table[h[n, k], {n, 0, 7}, {k, 0, n}] // TableForm

In the previous result, we notice  h(n, 0) - h(n, 1) = 1. We can verify the observation algebraically:
	\(TraditionalForm\`h(n, 0) - 
          h(n, 1)
 = \(D\_n - 
            n\ D\_\(n - 1\) = \(\(n!\)\ \(∑\+\(0 ≤ k ≤ n\)\((\(-1\))\)\^k\/\(k!\)\) - \(\(n(
                    n - 1)\)!\) \(∑\+\(0 ≤ k ≤ n - 1\)\((\(-1\))\)\^k\/\(k!\)\) = \((\(-1\))\)\^n\)\)\)

i.e., (Subscript[D, n]) = n (Subscript[D, n - 1]) + (((-1))^n)

Definition (Euler number)  
The number of permutations of the n objects having exactly k runs is denoted by ⟨({{n}, {k}})⟩.
	⟨({{n}, {k}})⟩=k⟨({{n}, {k}})⟩+(n + 1 - k)⟨({{n - 1}, {k - 1}})⟩, ⟨({{0}, {k}})⟩=I (k = 0) (i . e . , If[k == 0, 1, 0])

'crack'  breaks the given list into two parts at the specified position.

crack[l_, pos_]  := {Take[l, {1, pos}], 
        Take[l, {pos + 1, -1}]}

'cellDivision' explains the part, (n + 1 - k)⟨({{n - 1}, {k - 1}})⟩, of the recurrence relation for the Euler number.

\(cellDivision[l_, n_] := 
	
      Module[{dim, pos, id, s}, 
		dim\  = 
          Select[Transpose[{Range[Length[l]], 
						Length\  /@ \ 
                  l}], #[\([2]\)] =!= 1 &]; 
		pos = \(Flatten[#, 
                2] &\)[\(\(Outer[
                      List, {#1}, 
					Range[
                        2, #2]] &\) @@ # &\)\  /@ \ 
              dim]; 
		s = \(\(\(FlattenAt[#, 
                        id] &\)@
ReplacePart[l, 
                      crack[\(Insert[l, 
                            n, {id = #1, #2}]\)[\([#1]\)], #2], #] &\)
		 @@ # &\) /@ \ pos; 
				Prepend[s, 
          Prepend[l, {n}]]
							]\)

It is easy to see that the number of ways an element gets inserted at the end of each existing runs is k⟨({{n}, {k}})⟩ and at the positions other than the very end is (n + 1 - k)⟨({{n - 1}, {k - 1}})⟩. As an example, when an element x which is greater than 5 gets inserted into {{4,5},{2,3},{1}}, three different cases can occur, namely, {{x}, {4, 5}, {2, 3}, {1}},  {{4, x}, {5}, {2, 3}, {1}}, {{4, 5}, {2, x}, {3}, {1}}. However, {{4, 5}, {x}, {2, 3}, {1}}is accounted for the case where the element gets inserted at the end of each run.

l = {{4, 5}, {2, 3}, {1}};

cellDivision[l, x]

eulerRun[n_, k_] := {} /; n < k;
eulerRun[n_, 1] := {{Range[n]}};
eulerRun[n_, n_] := {Transpose[{Range[n, 1, -1]}]};
\(eulerRun[n_, k_] := 
	
      Join[
		Flatten[\(Table[
                Insert[#, n, {i, \(-1\)}], {i, 
                  k}] &\)\ 
			 /@ \ eulerRun[n - 1, k], 
          1], 
		Flatten[\(cellDivision[#, n] &\)\  /@ \ 
            eulerRun[n - 1, k - 1], 1]
	]\)

To obtain the set of permutations of {1, 2, 3, 4, 5} with 4 runs.

byRecursion = eulerRun[5, 4]

Length[%]

We verify that at least the number of runs are correct

Table[
      Length[eulerRun[n, k]] == Eulerian[n, k], {n, 6}, {k, n}]

To experimentally ascertain the correctness of the routine, 'eulerRun', we compare both the results from the routine and that from the brute force method of selecting the subsets with 4 runs from the total space.

totalSpace = Permutations[{1, 2, 3, 4, 5}];

allRuns = 
        Map[Split[#, #1 ≤ #2 &] &, totalSpace];

Select only the runs of size 4.

fromTotalSpace = Select[allRuns, Length[#] == 4 &]

Now change the form of the ouputs by using both methods from partition into the permutation

Flatten /@ byRecursion // Sort;

Flatten /@ fromTotalSpace // Sort;

% === %%