因为被一道毒瘤斯特林数题目卡掉,写一写斯特林数。

第一类斯特林数

        和组合数、排列数、卡特兰数等等相同,斯特林数是一个用在组合数学上的数,在很多计数问题上会用到。
        给出$n$个互异的数,问有多少中方案将它们组成$m$个圆排列?这个问题的答案就是第一类斯特林数,它的形式如下:

        这里的圆排列就是指指规定相对顺序不规定绝对顺序的排列,比如$ab$和$ba$是两个相同的圆排列,它们可以通过旋转相互得到。
        第一类斯特林数存在有符号数,这里只讨论无符号的那一种。对于无符号的第一类斯特林数,可以推导出它的递推公式。
        假若要求$S(n+1,m)$,可以考虑从$S(n,m)$和$S(n,m-1)$转移,对于$S(n,i),i<m-1$的情况,一定可以化到前两种情况中。从而可以得到递推公式:

        对于前者,单独的第$n+1$个元素可以独立形成一个圆排列,其数目就是$S(n,m)$。对于后者,将第$n+1$个元素插入到任意一个已有的$n$个元素左侧都可以得到一种方案,故为$nS(n,m-1)$。
        这个递推公式提供了$O(n^2)$求斯特林数的方法。
        显然,这样求的复杂度比较高,在$n$到$10^5$数量级时效率已经足以TLE,必须考虑其它的做法。
        第一类无符号斯特林数拥有生成函数形式:

        这个多项式展开后第$i$次方项系数就是$S(n,i)$。
        另外说一下,形如这样的多项式称为上升阶乘幂,即上升幂多项式,它常常用$x^{\overline n}$表示:

        现在问题就是如何展开这个多项式,考虑倍增,即求:

        问题转化为已知$F(x)$,求$F(x+k)$。假设$p_i$为$F(x)$的$i$次项系数,下面开始推式子:

        难搞的当然是右边那一个求和式,仔细一看有点像卷积,化一步,规定$G(x)=\frac {k^{-x}} {(-x)!}$,那么就是:

        显然是一个卷积,然后化成卷积熟悉的形式:

        考虑一下$a$和$b$的范围,为:

        $b$的次数为负数,这确实是很糟糕的。于是采用那个通用手段:多项式整体移动到非负范围。直接次数加上$n$即可:

        这样$b$的范围就是$i\leq b\leq n$,拿出来直接卷积即可。这里需要注意拿结果时需要再移回来。
        于是用$NTT$或者$FFT$在$O(nlogn)$下求出$F(x+k)$,再做一个多项式乘法,处理一下奇数时多余的项,这个工作就做完了。
        时间复杂度分析:

        复杂度就是$O(nlogn)$。模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include <bits/stdc++.h>

using namespace std;
#define N 300000
#define MOD 167772161
#define G 3
int tr[N << 2], fac[N], inv[N], n;
struct Pol {
int l, op[N << 2] = {0};
} pol1, pol2, pol3;

int qPow(int a, int b) {
int ans = 1, x = a;
while (b) {
if (b & 1)ans = 1ll * ans * x % MOD;
x = 1ll * x * x % MOD, b >>= 1;
}
return ans;
}


void NTT(int l, int *c, int type) {
for (int i = 0; i < l; i++)if (i < tr[i])swap(c[i], c[tr[i]]);
for (int mid = 1; mid < l; mid <<= 1) {
int wn = qPow(G, (MOD - 1) / (mid << 1));
if (type == -1)wn = qPow(wn, MOD - 2);
for (int len = mid << 1, j = 0; j < l; j += len) {
int w = 1;
for (int k = 0; k < mid; k++, w = 1ll * w * wn % MOD) {
int x = c[j + k], y = 1ll * w * c[j + mid + k] % MOD;
c[j + k] = (1ll * x + y) % MOD, c[j + mid + k] = (1ll * x - y) % MOD;
}
}
}
}

inline void multiple(Pol *ans, Pol *a, Pol *b) {
int l = 1, le = 0;
while (l <= a->l + b->l)l <<= 1, ++le;
for (int i = 0; i < l; i++) {
tr[i] = (tr[i >> 1] >> 1) | ((i & 1) << (le - 1));
if (i > a->l)a->op[i] = 0;
if (i > b->l)b->op[i] = 0;
}
NTT(l, a->op, 1), NTT(l, b->op, 1);
for (int i = 0; i < l; i++)a->op[i] = 1ll * a->op[i] * b->op[i] % MOD;
NTT(l, a->op, -1), l = qPow(l, MOD - 2), ans->l = a->l + b->l;
for (int i = 0; i <= ans->l; i++)ans->op[i] = 1ll * a->op[i] * l % MOD;
}

void solve(Pol *p, int x) {
if (x == 1) {
p->op[0] = 0, p->op[1] = 1, p->l = 1;
return;
}
int mid = x >> 1;
solve(p, mid);
for (int i = 0; i <= mid; i++)pol2.op[i] = 1ll * fac[i] * p->op[i] % MOD;
pol2.l = mid;
for (int i = 0; i <= mid; i++)pol3.op[i] = 1ll * qPow(mid, mid - i) * inv[mid - i] % MOD;
pol3.l = mid, multiple(&pol2, &pol2, &pol3);
for (int i = 0; i <= mid; i++)pol2.op[i] = 1ll * pol2.op[i + mid] * inv[i] % MOD;
pol2.l = mid, multiple(p, p, &pol2), pol3.op[0] = 0, pol2.op[x] = 0;
if (x % 2) {
for (int i = 0; i <= x - 1; i++)pol2.op[i] = 1ll * p->op[i] * (x - 1) % MOD;
for (int i = 1; i <= x; i++)pol3.op[i] = p->op[i - 1];
for (int i = 0; i <= x; i++)p->op[i] = (pol2.op[i] + pol3.op[i]) % MOD;
}
p->l = x;
}

int main() {
fac[0] = inv[0] = 1;
for (int i = 1; i <= 262145; i++)fac[i] = 1ll * i * fac[i - 1] % MOD, inv[i] = qPow(fac[i], MOD - 2);
cin >> n;
solve(&pol1, n);
for (int i = 0; i <= pol1.l; i++)cout << (pol1.op[i] + MOD) % MOD << " ";
return 0;
}

第二类斯特林数

        第二类斯特林数是关于集合划分的,$S(n,m)$表示将$n$个互异元素划分为$m$个非空集合的方案数,形式如下:

        第二类斯特林数拥有递推公式:

        前者是将多余的一个元素单独拿出来形成一个集合,后者则是放入任意一个集合中。这是容易理解的。
        递推式提供了$O(n^2)$求解的方法,下面是$O(nlogn)$的求法。
        第二类斯特林数的方案数公式:

        展开组合数,稍稍化一步:

        这是一个裸到不能再裸的卷积,直接卷一卷就可以在$O(nlogn)$下求出$S(n,0)->S(n,n)$。模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include <bits/stdc++.h>

using namespace std;
#define N 500000
#define MOD 167772161
#define G 3
int tr[N << 2], fac[N], inv[N], n;
struct Pol {
int l, op[N << 2] = {0};
} pol1, pol2;

int qPow(int a, int b) {
int ans = 1, x = a;
while (b) {
if (b & 1)ans = 1ll * ans * x % MOD;
x = 1ll * x * x % MOD, b >>= 1;
}
return ans;
}


void NTT(int l, int *c, int type) {
for (int i = 0; i < l; i++)if (i < tr[i])swap(c[i], c[tr[i]]);
for (int mid = 1; mid < l; mid <<= 1) {
int wn = qPow(G, (MOD - 1) / (mid << 1));
if (type == -1)wn = qPow(wn, MOD - 2);
for (int len = mid << 1, j = 0; j < l; j += len) {
int w = 1;
for (int k = 0; k < mid; k++, w = 1ll * w * wn % MOD) {
int x = c[j + k], y = 1ll * w * c[j + mid + k] % MOD;
c[j + k] = (1ll * x + y) % MOD, c[j + mid + k] = (1ll * x - y) % MOD;
}
}
}
}

inline void multiple(Pol *ans, Pol *a, Pol *b) {
int l = 1, le = 0;
while (l <= a->l + b->l)l <<= 1, ++le;
for (int i = 0; i < l; i++) {
tr[i] = (tr[i >> 1] >> 1) | ((i & 1) << (le - 1));
if (i > a->l)a->op[i] = 0;
if (i > b->l)b->op[i] = 0;
}
NTT(l, a->op, 1), NTT(l, b->op, 1);
for (int i = 0; i < l; i++)a->op[i] = 1ll * a->op[i] * b->op[i] % MOD;
NTT(l, a->op, -1), l = qPow(l, MOD - 2), ans->l = a->l + b->l;
for (int i = 0; i <= ans->l; i++)ans->op[i] = 1ll * a->op[i] * l % MOD;
}


int main() {
fac[0] = inv[0] = 1;
for (int i = 1; i <= 300000; i++)fac[i] = 1ll * i * fac[i - 1] % MOD, inv[i] = qPow(fac[i], MOD - 2);
cin >> n;
for (int i = 0; i <= n; i++)pol1.op[i] = (i % 2 == 0 ? 1 : -1) * inv[i];
for (int i = 0; i <= n; i++)pol2.op[i] = 1ll * qPow(i, n) * inv[i] % MOD;
pol1.l = pol2.l = n, multiple(&pol1, &pol1, &pol2);
for (int i = 0; i <= n; i++)cout << (1ll * pol1.op[i] + MOD) % MOD << " ";
return 0;
}

例题

        看一看那道斯特林数毒瘤题(其实也不是很毒瘤,只是我太菜):这里
        虽然最后也做出来了,不过是在比赛刚刚结束的时候:P
        如果我们能够求出以$AA..$或者$ABC..$开始的字符串数量,那么这个题就可以切掉了,求数量成了关键。
        假设现在前面是$ABA$(假设的),后面还有$n$个位待填,算一下$ABAA$开头的数量。首先后面需要进行集合划分,假如划分成$i$个,那么有$S(n,i)$(第二类斯特林数)种方案。
        由于有一些需要与前面已经划分好的结合,那么从中选出$p$个,这里$p\leq 2$且$p\leq n$,因为前面我们仅划分了两个集合,故$p\leq 2$。
        先组合再分配,于是方案数为:

        $k$是前面已经划分的数目,对于$ABA$为2。
        然后根据这个数跑一遍$DFS$就行了。注意$n=26$时爆long long,换用__int128。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#include <bits/stdc++.h>

using namespace std;

const int maxn = 28;
typedef __int128 ll;
ll Stirling[maxn][maxn];
ll C[30][30];
ll A[30][30];
int n;
ll k;

void init() {
Stirling[0][0] = 0;
Stirling[1][1] = 1;
for (ll i = 2; i < maxn; i++) {
for (ll j = 1; j <= i; j++) {
Stirling[i][j] = Stirling[i - 1][j - 1] + j * Stirling[i - 1][j];
}
}
C[1][0] = C[1][1] = 1;
for (int i = 2; i <= 26; i++) {
C[i][0] = C[i][i] = 1;
for (int j = 1; j < i; j++)C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
}
A[1][0] = A[1][1] = 1;
for (int i = 2; i <= 26; i++) {
A[i][0] = 1;
for (int j = 1, v = i; j <= i; j++, v--)A[i][j] = A[i][j - 1] * v;
}
}

ll get(int x, int y) {
if (x == 0)return 1;
ll p = 0;
for (int i = 1; i <= x; i++) {
ll s = 0;
for (int j = 0; j <= i && j <= y; j++)s += C[i][j] * A[y][j];
p += Stirling[x][i] * s;
}
return p;
}

void DFS(ll s, int x, int t) {
if (x == n + 1)return;
ll p;
for (int i = 1; i <= t; i++) {
p = get(n - x, t);
if (s <= p) {
printf("%c", 'A' + i - 1), DFS(s, x + 1, t);
return;
} else s -= p;
}
printf("%c", 'A' + t), DFS(s, x + 1, t + 1);
}

ll read() {
ll s = 0;
char e = getchar();
while (e < '-')e = getchar();
while (e > '-')s = (s << 1) + (s << 3) + (e & 15), e = getchar();
return s;
}

int main() {
init();
int t, p = 0;
scanf("%d", &t);
while (t--) {
scanf("%d", &n), k = read();
printf("Case #%d: ", ++p);
DFS(k, 1, 0);
printf("\n");
}
return 0;
}