Mathematica code for q-Stirling numbers

Question

In the paper [A new $q$-Analog of Stirling Numbers],(https://hal.archives-ouvertes.fr/hal-01372920/document)[PDF] J. Cigler defined $q$-Stirling numbers of the second kind as the following:

He considered the weight $w(\pi)$ for each partition $\pi$ of the set $A=\{0,\cdots,n-1\}$, and distinguished the part containing the zero element and called it $B_0$. Then,

$$w(\pi)=q^{\sum_{i\in B_0} i}$$ For each set $A$ of partitions, let $w(A)=\sum _{\pi\in A}w(\pi)$. Also, he considered $A_{n,k}$ to be the set of all partitions of $\{0,\cdots,n-1\}$ into $k$ non-empty parts. So $${n\brace k}=w(A_{n,k})$$

My question is:

What Mathematica code would generate these Stirling numbers of the second kind?

Answer 1

On second thought, it is much more convenient to use Theorem 4.5 of the paper as a definition, so

A[n_, k_, i_, j_] := Binomial[n, k + i] Binomial[n, k - j] - 
  Binomial[n, k + i + 1] Binomial[n, k - j - 1];

QStirlingS2Fast[n_, k_, q_] := 
  1/(1 - q)^(n - k) Sum[(-1)^i*
  A[n, k, i, j] q^Binomial[j + 1, 2] QBinomial[i, j, q], 
  {j, 0,k}, {i, j, n - k}];

This is of course much faster than the implementation via partitions

QStirlingS2[10, 5, 3] // RepeatedTiming
(* {2.00612, 3421737} *)

QStirlingS2Fast[10, 5, 3] // RepeatedTiming
(* {0.000111509, 3421737} *)

and works for large inputs as well

QStirlingS2Fast[1000, 20, 3] // N // RepeatedTiming  
(* {0.8409685, 1.465117*^8966} *)

Original post

Here is one possible implementation, which is a bit brute force. I am using Combinatorica` to obtain the set partitions, as well as the definition of the q-Stirling number of the second kind via the number of partition crossings $cr(\pi)$ given in the linked paper.

<< Combinatorica`

cr[pi_] := Cases[
  Flatten[
    Tuples /@ Map[
      Partition[#, 2, 1] &, 
      Subsets[Replace[pi, {_} :> Nothing, 1], {2}], 
    {2}], 
  1],
  {{i_, j_}, {k_, l_}} /; i < k < j < l || k < i < l < j
];


QStirlingS2[n_, k_, q_] := Sum[q^Length@cr[pi], {pi, KSetPartitions[Range@n, k]}]

This has the behavior for $q\to0$ given in eq.$(13)$,

Narayana[n_, k_] := 1/n*Binomial[n, k - 1] Binomial[n, k];
{ QStirlingS2[10, 4, q] /. q -> 0, Narayana[10, 4] }
(* {2520, 2520} *)

and also satisfies the identity of Theorem 4.5 in your reference

A[n_, k_, i_, j_] := Binomial[n, k + i] Binomial[n, k - j] - 
  Binomial[n, k + i + 1] Binomial[n, k - j - 1];

Theorem45[n_, k_, q_] := 
{
  QStirlingS2[n, k, q] (1 - q)^(n - k),
  Sum[(-1)^i*A[n, k, i, j] q^Binomial[j + 1, 2] QBinomial[i, j, q], 
  {j, 0, k}, {i, j, n - k}] 
}

Theorem45[6, 2, 3]
(* {1872, 1872} *)