Poly 杂记（三）

本章节以 Stirling Number 为中心展开。

小式子

\(p\) 是素数时

\[x^{\overline p} = x^p - x \pmod p \]

换句话说，由于

\[x^{\overline p} = \sum_{i=0}^p {p \brack i} x^i \]

所以在模 \(p\) 意义下

\[{p \brack p} = 1,{p \brack 1} = -1,{p \brack 2}= {p \brack 3} = \cdots = {p \brack p-1} =0 \]

证明：

根据 威尔逊定理，我们可以知道 \({p \brack 1} = (p-1)! \equiv -1 \pmod p\)，前两者容易证明。

其余部分我们考虑直接构造两个多项式

\[\begin{aligned} S(x) = & x(x+1)(x+2) \cdots (x+p-1) = x^{\overline p} \\ F(x) = & x^p - x \end{aligned} \]

容易发现两者都是以 \(0 \sim p-1\) 为根的，于是我们直接做差可以得到

\[G(x) = S(x) - F(x) \]

有 \(p\) 个 \(0\) 点，而 \(G(x)\) 次数 \(\lt p\)（最高次项 \(x^p\) 抵消了）。

所以 \(G(x) = 0\)（根据拉格朗日定理/多项式推理法）。得证。

这个式子在 QOJ 1285 中有应用。

另外一个神秘之师

\[k ^{\underline i} \binom n k = n^{\underline i} \binom{n-i} {k-i} \]

直接展开即可得证。

斯特林数的四种求法

有一回对我说道，“你读过书么？”我略略点一点头。他说，“读过书，……我便考你一考。组合数学里的斯特林数，怎样说的？”我想，AKIOI的人，我也配答么？便回过脸去，不再理会。田乙己等了许久，很恳切的说道，“不会罢？……我教给你，记着！这些数应该记着。将来做OIer的时候，省选要用。”我暗想我和省选的等级还很远呢，而且我们考纲里也没有斯特林数；又好笑，又不耐烦，懒懒的答他道，“谁要你教，不是将n个元素划分成m个集合或环排列的方案数嘛？”田乙己显出极高兴的样子，将两个指头的长指甲敲着鼠标，点头说，“对呀对呀！……斯特林数有四样求法，你知道么？”我愈不耐烦了，努着嘴走远。田乙己刚打开typora，想在电脑上打公式，见我毫不热心，便又叹一口气，显出极惋惜的样子。

第一类斯特林数 - 行

求 \({n \brack 0} \sim {n \brack n}\)，转求 \(x^{\overline n}\) 的系数。

定义 \(F_n(x) = x^{\overline n}\)，则

\[F_{2n}(x) = F_n(x) F_n(x+n) \]

于是多项式平移 + 倍增即可。时间复杂度 \(\mathcal O(n \log n)\)。

模板：P5408 第一类斯特林数·行。代码。

Code

poly a,b;

poly Move(poly F,int c){
  int n=(int)F.size();
  a.resize(n),b.resize(n);
  for(int i=0;i<n;i++) a[i]=mul(F[n-i-1],fac[n-i-1]),b[i]=mul(ifac[i],qpow(c,i));
  a=a*b;
  for(int i=0;i<n;i++) F[i]=mul(ifac[i],a[n-1-i]);
  return F;
}

map<int,poly> mp;

poly solve(int n){
  if(mp.count(n)) return mp[n];
  if(n==1) return {0,1};
  return mp[n]=solve(n>>1)*Move(solve(n-(n>>1)),n>>1);
}

int main(){
  /*2024.7.22 H_W_Y P5408 第一类斯特林数·行 多项式平移*/
  cin>>n;init(n);
  if(n==0) cout<<1,exit(0);
  F=solve(n);
  for(int i=0;i<=n;i++) cout<<F[i]<<' ';
  return 0;
}

第二类斯特林数 - 行

求 \({n \brace 0} \cdots {n \brace n}\)。

我们有如下恒等式

\[m^n = \sum_{i=0}^m \binom m i i! {n \brace i} \]

这个式子也是简单的，就是考虑看成把 \(n\) 个互不相同的球放入 \(m\) 个互不相同的盒子里面，后者就是枚举那些盒子是空的。

于是我们对上面的式子进行二项式反演

\[\begin{aligned} m! {n \brace m} & = \sum_{i=0}^m \binom m i (-1)^{m-i} i^n \\ {n \brace m} & = \sum_{i=0}^m \frac {\binom m i(-1)^{m-i} i^n} {m!} \\ & =\sum_{i=0}^m \frac {i^n} {i!} \frac {(-1)^{m-i}} {(m-i)!} \end{aligned} \]

最后的式子是一个卷积形式，所以我们可以用多项式乘法完成 第二类斯特林数行 的求法。代码。

Code

void Mul(ull *F,ull *G,int n,int m){
  int P=1,I;
  for(;P<n+m;P<<=1);
  Init(P),I=qpow(P);
  NTT(F,P,1),NTT(G,P,1),Cro(F,G,P),NTT(F,P,-1);
  for(int i=0;i<P;i++) F[i]=mul(F[i],I);
}

int main(){
  /*2024.6.27 H_W_Y P5395 第二类斯特林数·行*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n;init();

  for(int i=0;i<=n;i++){
	F[i]=mul(qpow(i,n),ifac[i]);
    if(i&1) G[i]=dec(0,ifac[i]);
    else G[i]=ifac[i];
  }
  Mul(F,G,n+1,n+1);
  for(int i=0;i<=n;i++) cout<<F[i]<<' ';
  return 0;
}

我们还有另外一种求法，即将 \(x^n\) 展开成下降幂形式

\[x^n = {n \brace 0} + {n \brace 1}x ^{\underline 1} + {n \brace 2} x^{\underline 2} + \cdots \]

转化为求 \(0 \sim n\) 的点值，再用之前下降幂系数的 Trick \(\mathcal O(n \log n)\) 求得。（在杂记（二）里面）。

第一类斯特林数 - 列

接下来我们考虑用生成函数来求。

求 \({0 \brack k} \sim {n \brack k}\)。

设 \(n\) 个元素放入一个圆排列，也就是 \({n \brack 1} = (n-1)!\)，

于是我们有它的 EGF 是 \(\hat A(x) = \sum_{n \gt 0} \frac {x^n} n\)。

那么放入 \(k\) 个圆排列的 EGF 就是

\[\frac 1 {k!} \hat A(x) \]

前面乘上 \(\frac 1 {k!}\) 是因为圆排列是无序的。

直接用多项式快速幂解决即可。时间复杂度 \(\mathcal O(n \log n)\)。

模板：P5409 第一类斯特林数·列。代码。

Code

poly Qpow_plus(poly F,int k){
  int n=(int)F.size(),id=0;
  poly A;

  for(int i=0;i<n;i++) if(F[i]!=0){
	int I=qpow(F[i]);
	A.resize(n-id);
	for(int j=i;j<n;j++)
	  A[j-i]=mul(F[j],I);
	id=i;break;
  }
  A=Qpow(A,k);

  int nw=qpow(F[id],k);
  id=mul(id,k);
  for(int i=0;i<n;i++) F[i]=0;
  for(int i=id;i<n&&i<(int)A.size()+id;i++) F[i]=mul(A[i-id],nw);
  return F;
}

int main(){
  /*2024.7.22 H_W_Y P5409 第一类斯特林数·列*/
  cin>>n>>m;init(N-1);
  F.resize(n+1);
  for(int i=1;i<=n;i++) F[i]=inv[i];
  F=Qpow_plus(F,m);
  for(int i=0;i<=n;i++) cout<<mul(F[i],mul(fac[i],ifac[m]))<<' ';
  return 0;
}

第二类斯特林数 - 列

求 \({0 \brace k} \sim {n \brace k}\)。

还是考虑生成函数。

同样的，我们把 \(n\) 个元素放入一个集合中，那么 \({n \brace 1} = 1\)，于是其 EGF 为

\[\sum_{n \gt 0} \frac {x^n} {n!}= e^x-1 \]

那么放入 \(k\) 个集合中的答案就是

\[\frac 1 {k!} (e^x-1)^k \]

于是就做完了，时间复杂度 \(\mathcal O(n \log n)\)。

模板：P5396 第二类斯特林数·列。代码。

Code

int main(){
  /*2024.7.22 H_W_Y P5396 第二类斯特林数·列*/
  cin>>n>>m;init(N-1);
  F.resize(n+1);
  for(int i=1;i<=n;i++) F[i]=ifac[i];
  F=Qpow_plus(F,m);
  for(int i=0;i<=n;i++) cout<<mul(F[i],mul(fac[i],ifac[m]))<<' ';
  return 0;
}

和上一道题的代码 \(\Delta = 3\)，还是比较牛的。

P4091 [HEOI2016/TJOI2016] 求和 - 组合意义

现在他想计算这样一个函数的值:

\[f(n)=\sum_{i=0}^n\sum_{j=0}^i S(i,j)\times 2^j \times (j!) \]

S(i, j)表示第二类斯特林数，递推公式为:

\(S(i, j) = j \times S(i - 1, j) + S(i - 1, j - 1), 1 \le j \le i - 1\)。

边界条件为：\(S(i, i) = 1(0 \le i), S(i, 0) = 0(1 \le i)\)

你能帮帮他吗?

\(1 \leq n \leq 10^5\)。

我们现在考虑求解

\[g_n =\sum_{i=0}^n {n \brace i} 2^i i! \]

考虑其组合意义：

把 \(n\) 个球放入 \(i\) 个 box 中的方案数，每个 box 不相同且对 box 进行黑白染色。

上面的式子枚举的是球，我们考虑枚举 box，即枚举最后一个 box 放了多少个球，于是有

\[g_n = \sum_{i=1}^n 2 \binom n i g_{n-i} \]

用 EGF 的常用套路展开一下就可以得到

\[\frac {g_n} {n!} = \sum_{i=1}^n \frac 2 {i!} \frac {g_{n-i}} {(n-i)!} \]

于是设

\[G(x) = \sum \frac {g_n} {n!},F(x) = \sum_{n \gt 0} \frac 2 {n!} \]

上式就是

\[\begin{aligned} G(x) = & G(x) F(x)+1 \\ \therefore G(x) = & \frac 1 {1-F(x)} \end{aligned} \]

直接多项式求逆即可，时间复杂度 \(\mathcal O(n \log n)\)。代码。

Code

int main(){
  /*2024.7.23 H_W_Y P4091 [HEOI2016/TJOI2016] 求和 Striling Number*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n;F.resize(n+1),init(n);
  for(int i=1;i<=n;i++) F[i]=dec(0,mul(2,ifac[i]));
  F[0]=1;Inv(F);
  for(int i=0;i<=n;i++) add(ans,mul(F[i],fac[i]));
  cout<<ans;
  return 0;
}

P6620 [省选联考 2020 A 卷] 组合数问题 - 下降幂/Snake Oil

众所周知，小葱同学擅长计算，尤其擅长计算组合数。小葱现在希望你计算

\[\left(\sum_{k=0}^{n}f(k)\times x^k\times \binom{n}{k}\right)\bmod p \]

的值。其中 \(n\), \(x\), \(p\) 为给定的整数，\(f(k)\) 为给定的一个 \(m\) 次多项式 \(f(k) = a_0 + a_1k + a_2k^2 + \cdots + a_mk^m\)。\(\binom{n}{k}\) 为组合数，其值为 \(\binom{n}{k} = \frac{n!}{k!(n-k)!}\)。

\(1\le n, x, p \le 10^9, 0\le a_i\le 10^9, 0\le m \le \min(n,1000)\)。

由于涉及求和，所以我们考虑 下降幂，即把 \(f\) 转化成下降幂形式。

考虑 \(f(k) = k ^{\underline m}\) 的情况，那么其它的情况我们都可以解决了。

于是式子就是

\[\sum_{k=0}^n k ^{\underline m} x^k \binom n k \]

考虑一个神秘式子

\[k ^{\underline i} \binom n k = n^{\underline i} \binom{n-i} {k-i} \]

进而我们可以优化上面的式子了

\[\begin{aligned} & \sum_{k=0}^n k ^{\underline m} x^k \binom n k \\ = & \sum_{k=0}^n x^k n^{\underline m} \binom {n-m} {k-m} \\ = & n^{\underline m} \sum_{k=0}^{n-m} x^{m+k} \binom {n-m} k \\ = & n^{\underline m} x^m \sum_{k=0}^{n-m} x^k \binom {n-m} k \\ = & n^{\underline m} x^m (x+1)^{n-m} \end{aligned} \]

如果我们把 \(f(x)\) 表示成

\[f(x) = \sum_{i=0}^m b_i x^{\underline i} \]

那么可以得到答案就是

\[ans= \sum_{i=0}^m b_i n^{\underline i} x^i (x+1)^{n-i} \]

由于 \(m \le 1000\)，所以我们可以用 \(\mathcal O(n^2)\) 的时间得到 \(f\) 的下降幂形式。

但是由于 \(p\) 是任意的整数，所以我们不能涉及求逆操作，于是普通幂到下降幂的转化可以直接用 第二类斯特林数 完成，也就是

\[\begin{aligned} & \sum_{i=0}^m a_i x^i \\ = & \sum_{i=0}^m a_i \sum_{j=0}^i {i \brace j} x^{\underline j} \\ = & \sum_{i=0}^m x^{\underline i} \sum_{j=i}^m {j \brace i} a_j \end{aligned} \]

于是用递推式预处理出第二类斯特林数就可以完成了。

时间复杂度 \(\mathcal O(m^2)\)。代码。

Code

void init(int n){
  S[0][0]=1;
  for(int i=1;i<=n;i++)
    for(int j=1;j<=i;j++)
	  S[i][j]=adc(mul(j,S[i-1][j]),S[i-1][j-1]);
}

int main(){
  /*2024.7.23 H_W_Y P6620 [省选联考 2020 A 卷] 组合数问题*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>x>>H>>m;init(m);
  for(int i=0;i<=m;i++) cin>>a[i];
  for(int i=0;i<=m;i++) for(int j=i;j<=m;j++) add(b[i],mul(S[j][i],a[j]));
  for(int i=0,p1=1,p2=1;i<=m;i++){
	add(ans,mul(mul(b[i],p1),mul(p2,qpow(x+1,n-i))));
	p1=mul(p1,n-i),p2=mul(p2,x);
  }
  cout<<ans;
  return 0;
}

不用 下降幂？

这道题我们还可以用另外一种方法推一推，似乎叫 Snake Oil。

直接把式子展开

\[\begin{aligned} & \sum_{k=0}^n \binom n k x^k \sum_{i=0}^m a_i k^i \\ = & \sum_{i=0}^m a_i \sum_{k=0}^n \binom n k x^k k ^i \end{aligned} \]

我们设

\[f_i = \sum_{k=0}^n \binom n k x^k k^i \]

于是需要敏锐地观察到 \(f\) 的 EGF 是可求的：

\[\begin{aligned} f(z) = & \sum_{i=0}^{\infty} \frac {z^i} {i!} f_i \\ = & \sum_{i=0}^{\infty} \frac{z^i} {i!} \sum_{k=0}^n \binom n k x^k k^i \\ = & \sum_{k=0}^n \binom n k x^k \sum_{i=0}^{\infty} \frac {z^i k^i} {i!} \\ = & \sum_{k=0}^n \binom n k (xe^z) ^k \\ = & (1+xe^z)^n \end{aligned} \]

于是这个可以用多项式快速幂解决。但是实际上在这道题中不适用。

P6667 [清华集训2016] 如何优雅地求和 - 下降幂/PGF

有一个多项式函数 \(f(x)\)，最高次幂为 \(x^m\)，定义变换 \(Q\)：

\[Q(f,n,x)=∑_{k=0}^{n}f(k)\binom nkx^k(1−x)^{n−k} \]

现在给定函数 \(f\) 和 \(n,x\)，求 \(Q(f,n,x)\bmod998244353\)。

出于某种原因，函数 \(f\) 由点值形式给出，即给定 \(a_0,a_1,⋯,a_m\) 共 \(m+1\) 个数，\(f(x)=a_x\)。可以证明该函数唯一。

\(1≤n≤10^9\)，\(1≤m≤2×10^4\)，\(0≤a_i,x<998,244,353\)。

出于某种原因？给连续点值该是提醒你用下降幂了/qd

和上一题基本一样，我们还是先考虑 \(f(k) = k^{\underline m}\) 的情况，这道题只是多了一个 \((1-x)^{n-k}\)，那么

\[\begin{aligned} & \sum_{k=0}^n k^{\underline m} \binom n k x^k (1-x)^{n-k} \\ = & \sum_{k=0}^n n^{\underline m} \binom {n-m} {k-m} x^k (1-x)^{n-k} \\ = & n^{\underline m} \sum_{k=0}^{n-m} \binom {n-m} k x^{m+k} (1-x)^{n-m-k} \\ = & n^{\underline m} x^m \sum_{k=0}^{n-m} \binom {n-m} k x^k (1-x)^{n-m-k} \\ = & n^{\underline m} x^m (x+1-x)^{n-m} \\ = & n^{\underline m} x^m \end{aligned} \]

非常优美的式子啊，于是答案就是

\[\sum_{i=0}^m b_i n^{\underline i} x^i \]

现在问题就变成了求下降幂系数 \(b_i\) 了，这是杂记（二）中的经典问题，于是直接 \(\mathcal O(n \log n)\) 卷积就做完了。代码。

Code

int main(){
  /*2024.7.23 H_W_Y P6667 [清华集训2016] 如何优雅地求和 下降幂*/
  cin>>n>>m>>x;init(m);

  F.resize(m+1),G.resize(m+1);
  for(int i=0;i<=m;i++){
	cin>>F[i],F[i]=mul(F[i],ifac[i]);
    if(i&1) G[i]=dec(0,ifac[i]);
	else G[i]=ifac[i];
  }
  F=F*G;

  for(int i=0,p1=1,p2=1;i<=m;i++)
	add(ans,mul(F[i],mul(p1,p2))),p1=mul(p1,n-i),p2=mul(p2,x);

  cout<<ans;
  return 0;
}

这道题还有一种用 概率生成函数 PGF 的方法。

我们稍微改一下原式的变量名 \(x \to p\)，于是可以得到

\[\sum_{k=0}^n f(k) \binom n k p^k (1-p)^{n-k} \]

这个不是和一个事件发生的概率看起来这么像？概率分布？

考虑一种实验，有 \(p\) 的概率成功，有 \(1-p\) 的概率失败。

而 \(X\) 是一种随机变量，表示 \(n\) 次实验中成功的次数，于是我们有

\[P_r(X = k) = \binom n k p^k (1-p)^{n-k} \]

于是原式就是

\[Q = \sum_{k=0}^n f(k) P_r(X=k) = E(f(X)) \]

同样的，我们设 \(f(k) = x^{\underline m}\)，于是我们可以得到 \(E(x^{\underline m})\) 的 概率生成函数 PGF。

也就是一次实验的 PGF 是 \((pz + (1-p))\)，则 \(X\) 的 PGF 是 \((pz+(1-p))^n\)。

然后你需要观察一个结论：\(X\) 的 PGF 为 \(G(x)\) 则

\[E(x^{\underline m}) = G^{(m)}(1) = n ^{\underline m } p^m \]

证明就是发现

\[G^{(k)}(x) = g_k k^{\underline k} + {g_{k+1}} (k+1)^{\underline k} + \cdots \]

所以我们有

\[G^{(k)} (1) = \sum_{i \ge k} g_i i^{\underline k} = \sum_{i \ge k} i^{\underline k} P_r(x=i) = E(x^{\underline k}) \]

证毕。而下一步的化简更为显然。

于是你又可以像上一种方法那样做了。

感觉这种方法不是很优？你需要观察到很多东西啊/jk

但是，这里有一道几乎一模一样的例题：CF1278F Cards。

直接把 \(x^k\) 用 Stirling 数展开成下降幂形式，然后用上述方法求解就可以了，这里 \(p=\frac 1 m\)。代码。

Code

#include <bits/stdc++.h>
using namespace std;

const int N=5e3+5,H=998244353;
int n,m,k,p,s[N][N],ans;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

void init(){
  s[0][0]=1;
  for(int i=1;i<=k;i++)
    for(int j=1;j<=i;j++) s[i][j]=adc(s[i-1][j-1],mul(s[i-1][j],j));
  p=qpow(m);
}

int main(){
  /*2024.8.1 H_W_Y CF1278F Cards PGF*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>m>>k;init();
  for(int i=0,p1=1,p2=1;i<=k;i++)
    add(ans,mul(s[k][i],mul(p1,p2))),p1=mul(p1,n-i),p2=mul(p2,p);
  cout<<ans;
  return 0;
}

QOJ 1285 Stirling Number

给定四个整数 \(n,l,r,p\)，求

\[\left( \sum_{k=l}^r {n \brack k} \right) \bmod p \]

\(1 \le n \le 10^{18},0 \le l \le r \le n,2 \le p \le 10^6\)，保证 \(p\) 是质数。

好题！

由于一行的第一类斯特林数是由 \(x^{\overline n}\)，所以我们考虑优化这个部分，设

\[S_n(x)= x^{\overline n} \]

发现我们需要用 \(\bmod p\) 的性质，由于我们知道

\[x^{\overline p} = x^p -1 \pmod p \]

所以设 \(n=qp+r\)，于是有

\[\begin{aligned} S_n(x) = & x^{\overline n} \\ = & \prod_{i=0}^{n-1} (x+i) \\ = & \prod_{i=0}^{q-1} (x+ip)(x+ip+1)(x+ip+2)\cdots (x+ip+p-1)\prod_{i=0}^{r} (x+qp+i) \\ = & \prod_{i=0}^{q-1} x^{\overline p} \prod_{i=0}^{r} (x+i) \\ = & (x^p-x)^q S_r(x) \\ = & x^q (x^{p-1}-1)^q S_r(x) \end{aligned} \]

那么我们考虑 \(n \brack k\) 的表示

\[\begin{aligned} {n \brack k} = & [x^k] S_n(x) \\ = & [x^k] x^q (x^{p-1}-1)^q S_r(x) \\ = & [x^{k-q}] \prod_{i=0}^{q} \binom q i x^{(p-1)i} (-1)^{q-i} \sum_{j=1}^r {r \brack j} x^j \end{aligned} \]

注意，这里我们钦定 \(r \neq 0\)，也就是 \(r=0\) 的情况在后面特殊处理，于是枚举的 \(i\) 是从 \(1\) 开始的。

发现系数一定是形如 \((p-1)i + j\) 的形式，于是我们设 \(k-q= a(p-1)+b\)，其中 \(1 \le b \le r\)，于是我们有

\[{n \brack k} = \binom q a (-1)^{q-a} {r \brack b} \]

于是把最初的式子差分之后，我们的目标是计算第一类斯特林数的前缀和，也就是可以得到

\[\sum_{i=1}^M {n \brack i} = \sum_{i=1}^M \binom q a (-1)^{q-a} {r \brack b} \]

由于 \(p\) 很小，所以我们考虑枚举 \(b\)，设 \(M-q = m (p-1) + k(1 \le k \le p-1)\)，于是可以得到

\[\begin{aligned} \sum_{i=1}^M {n \brack i} = & \sum_{i=1}^k {r \brack i} (-1)^q \sum_{j=0}^m \binom q j (-1)^j + \sum_{i=k+1}^r {r \brack i} (-1)^q \sum_{j=0}^{m-1} \binom q j (-1)^j \end{aligned} \]

对于

\[\sum_{i=0}^m \binom q i (-1)^i \]

我们可以直接用上指标反转之后平行求和，也就是可以得到

\[\begin{aligned} \sum_{i=0}^m \binom q i (-1)^i = & \sum_{i=0}^m \binom {i-q-1} i = \binom {m-q} m = (-1)^m \binom {q-1}{m} \end{aligned} \]

那么原式就是

\[\begin{aligned} \sum_{i=1}^M {n \brack i} = & (-1)^q \left( \sum_{i=1}^{k} {r \brack i} (-1)^m \binom {q-1} m + \sum_{i=k+1}^r {r \brack i} (-1)^{m-1} \binom {q-1} {m-1} \right) \\ = & (-1)^{q+m} \left( \sum_{i=0}^k {r \brack i} \binom {q-1} m - \sum_{i=k+1}^r {r \brack i} \binom {q-1} {m-1} \right) \\ = & (-1)^{q+m} \left( \sum_{i=0}^k {r \brack i} \left( \binom {q-1} m + \binom {q-1} {m-1} \right) - \sum_{i=0}^r {r \brack i} \binom {q-1} {m-1} \right) \\ = & (-1)^{q+m} \left( \binom {q} m \sum_{i=0}^k{r\brack i} - \binom {q-1} {m-1} r! \right) \end{aligned} \]

那么现在问题又变成了求

\[\sum_{i=0}^k {r \brack i} \]

一个最简单的做法是直接把这一行求出来，也就是用 任意模数 NTT 直接做，时间复杂度 \(\mathcal O(p \log p)\)，但是我们存在 \(\mathcal O(p)\) 的线性做法，而且较为好写。

其中 \(r \lt p\)，根据 EI 的这个（小质数单点S1） 我们可以用 \(f(1),f(g),f(g^2),\cdots,f(g^{p-2})\) 去 IDFT 出系数，这里的 \(f(x) = x^{\overline r}\)，其中 \(g\) 是 \(p\) 意义下的原根。

也就是说，我们有式子

\[{r \brack i} = \frac 1 {p-1} \sum_{j=0}^{p-2} f(g^j) g^{-ij} \]

于是问题的式子就是

\[\begin{aligned} \sum_{i=0}^k {r \brack i} = & \sum_{i=0}^k \frac {1} {p-1} \sum_{j=0}^{p-2} f(g^j) g^{-ij} \\ = & - \sum_{i=0}^k \sum_{j=1}^{p-2} f(g^j) g^{-ij} - (k+1) f(1) \\ = & - \sum_{j=1}^{p-2} f(g^j) \sum_{i=0}^k g^{-ij} - (k+1) f(1) \\ = & - \sum_{j=1}^{p-2} f(g^j) \frac {g^{-(k+1)j} -1} {g^{-j} -1} - (k+1)f(1) \end{aligned} \]

就可以在 \(\mathcal O(p)\) 的线性时间复杂度之内完成了。

我们似乎还落下了一个 \(r=0\) 的情况，也来小推一手，发现是容易的

\[\begin{aligned} {n \brack k} = & [x^k] S_n(x) \\ = & [x^k] x^q (x^{p-1}-1)^q S_r(x) \\ = & [x^{k-q}] \prod_{i=0}^{q} \binom q i x^{(p-1)i} (-1)^{q-i} \end{aligned} \]

也就是说这里需要 \((p-1) \mid k-q\)，同样设 \(M-q = m(p-1) +k\)，那么答案就是

\[\begin{aligned} \sum_{i=0}^M {n \brack i} = & \sum_{i=0}^m \binom q i (-1)^{q-i} \\ = & (-1)^{q+m} \binom {q-1} m \end{aligned} \]

于是这道题就做完了，时间复杂度 \(\mathcal O(p)\)。代码。

Code

#include <bits/stdc++.h>
using namespace std;
#define int long long

const int N=2e6+5;
int L,R,n,q,p,r,m,k,fac[N],ifac[N],f[N],prm[N],ct=0,g,pw[N],ipw[N];

int adc(int a,int b){return a+b>=p?a+b-p:a+b;}
int dec(int a,int b){return a<b?a-b+p:a-b;}
int mul(int a,int b){return 1ll*a*b%p;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=p-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

bool chk(int x){
  if(qpow(x,p-1)!=1) return 0;
  for(int i=1;i<=ct;i++)
    if(qpow(x,(p-1)/prm[i])==1) return 0;
  return 1;
}

int findrt(){
  for(int i=2;i<p;i++) if(chk(i)) return i;
  return -1;
}

void init(){
  fac[0]=1;
  for(int i=1;i<p;i++) fac[i]=mul(fac[i-1],i);
  ifac[p-1]=qpow(fac[p-1]);
  for(int i=p-1;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
  
  int x=p-1;
  for(int i=2;i*i<=x;i++){
	if(x%i==0){
	  prm[++ct]=i;
	  while(x%i==0) x/=i;
	}
  }
  if(x!=1) prm[++ct]=x;

  g=findrt();

  pw[0]=1,f[0]=fac[r];
  for(int i=1;i<=p-2;i++)
	pw[i]=mul(pw[i-1],g),f[i]=mul(fac[pw[i]+r-1],ifac[pw[i]-1]);
}

int C(int n,int m){
  if(m<0||n<m) return 0;
  return mul(fac[n],mul(ifac[m],ifac[n-m]));
}

int binom(int n,int m){
  if(m<0||n<m) return 0;
  if(m==0) return 1;
  return mul(C(n%p,m%p),binom(n/p,m/p));
}

int calc(int k){
  if(k>=r) return fac[r];
  int res=dec(0,mul(k+1,f[0]));
  for(int j=1;j<=p-2;j++)
    del(res,mul(f[j],mul(dec(pw[(p-1-j)*(k+1)%(p-1)],1),qpow(dec(pw[p-1-j],1)))));
  return res;
}

int SOLVE(int M){
  if(!M||M<q) return 0;
  m=(M-q)/(p-1),k=(M-q)%(p-1);

  if(!r){
	if((q+m)&1) return dec(0,binom(q-1,m));
	return binom(q-1,m);
  }

  if(!k) k=p-1,--m;
  int res=0;

  res=dec(mul(binom(q,m),calc(k)),mul(binom(q-1,m-1),fac[r]));
  if((q+m)&1) return dec(0,res);
  return res;
}

signed main(){
  /*2024.7.25 H_W_Y QOJ 1285 Stirling Number*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>L>>R>>p;
  q=n/p,r=n%p;init();
  cout<<dec(SOLVE(R),SOLVE(L-1));
  return 0;
}

CF932E Team Work - Stirling

给定 $ n,k $，求：

\[\sum_{i=1}^n\binom n i \times i^k \]

$ 1 \leq k \leq 5000,1 \leq n \leq 10^9 $

求和？直接展开成下降幂形式，即用第二类斯特林数展开

\[\begin{aligned} & \sum_{i=1}^n \binom n i \sum_{j=0}^k {k \brace j} i^{\underline j} \\ = & \sum_{j=0}^k {k \brace j} \sum_{i=0}^n \binom n i i^{\underline j} \\ = & \sum_{j=0}^k {k \brace j} \sum_{i=0}^n \binom {n-j} {i-j} n ^{\underline j} \\ = & \sum_{j=0}^k {k \brace j} n^{\underline j} \sum_{i=0}^{n-j} \binom {n-j} i \\ = & \sum_{j=0}^k {k \brace j} n^{\underline j} 2^{n-j} \end{aligned} \]

于是直接 \(k^2\) 求第二类斯特林数然后暴力求就可以了。

时间复杂度 \(\mathcal O(k^2)\)。代码。

Code

#include <bits/stdc++.h>
using namespace std;

const int N=5e3+5,H=1e9+7;
int s[N][N],n,k,inv2,ans;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

void init(){
  s[0][0]=1;
  for(int i=1;i<=k;i++)
    for(int j=1;j<=i;j++) s[i][j]=adc(s[i-1][j-1],mul(s[i-1][j],j));
  inv2=qpow(2);
}

int main(){
  /*2024.8.1 H_W_Y CF932E Team work Stirling*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>k;init();
  for(int i=0,p1=qpow(2,n),p2=1;i<=k;i++) add(ans,mul(s[k][i],mul(p1,p2))),p1=mul(p1,inv2),p2=mul(p2,n-i);
  cout<<ans;
  return 0;
}

CF717A Festival Organization - Fibonacci + Stirling + 扩域

一个合法的串定义为：长度在 \([l,r]\) 之间，且只含 0,1，并且不存在连续 \(2\) 个或更多的 \(0\)。

现在要选出 \(k\) 个长度相同的合法的串，问有几种选法，答案模 \(10^9+7\)。

不是所有 \(a^b\) 都有欧拉定理。好题。

首先，考虑普通 dp，设 \(f_{i,0/1}\) 表示长度为 \(i\) 最后一位是 \(0\) 还是 \(1\) 的方案数，那么容易有一下转移式子

\[f_{i,0} = f_{i-1,1} ,f_{i,1} = f_{i-1,0} + f_{i-1,1} \]

一看就有斐波那契数列的感觉，进而，我们发现

\[f_{i,0}+f_{i,1} = F_{i+2} \]

其中 \(F_i\) 表示 Fibonacci 数列的第 \(i\) 项。

那么其实答案就是

\[ans = \sum_{i=l+2}^{r+2} \binom {F_{i}} k \]

我们令 \(L=l+2,R=r+2\)，把组合数转化成下降幂形式，就可以得到

\[ans = \frac1 {k!} \sum_{i=L}^R {F_i}^{\underline k} \]

在这里，下降幂明显是不好处理的，所以我们考虑用 带符号第一类斯特林数 转化成幂形式

\[\begin{aligned} ans = & \frac 1 {k!} \sum_{i=L}^R {F_i}^{\underline k} \\ = & \frac 1 {k!} \sum_{i=L}^R \sum_{j=0}^k (-1)^{k-i} {k \brack j} {F_i}^k \end{aligned} \]

对于 Fibonacci，我们能说什么？

容易知道其生成函数为 \(\frac x {1-x-x^2}\)，那么我们有其 通项公式

\[F_n = \frac {x^n - y^n} {\sqrt 5},x=\frac {1+\sqrt 5} 2,y=\frac {1- \sqrt 5} 2 \]

这里我们令 \(a = \frac 1 {\sqrt 5},b = - \frac 1 {\sqrt 5}\)，那么其实

\[F_n = ax^n + by^n \]

把它带入上面的式子，可以得到

\[\begin{aligned} ans = & \frac 1 {k!} \sum_{i=0}^k (-1)^{k-i} {k \brack i} \sum_{j=L}^R \left( ax^j + by^j \right)^i \\ = & \frac 1 {k!} \sum_{i=0}^k (-1)^{k-i} {k \brack i } \sum_{j=L}^R \sum_{p=0}^i \binom i p a^p (x^j)^p b^{1-p} (y^j)^{1-p} \\ = & \frac 1 {k!} \sum_{i=0}^k (-1)^{k-i} {k \brack i} \sum_{p=0}^i \binom i p a^p b^{1-p} \sum_{j=L}^R \left( x^py^{1-p} \right)^j \end{aligned} \]

最后的那一项可以用等比数列求和公式完成。于是这样我们似乎就做完了。

真的做完了吗？发现并非如此，因为 \(5\) 在 \(10^9+7\) 的意义下没有二次剩余。

所以我们考虑 扩域，把每个数都表示成 \(a+b \sqrt 5\)，这样运算法则和复数类似，就可以完成了。

时间复杂度 \(\mathcal O(n^2)\)。代码。

Code

#include <bits/stdc++.h>
using namespace std;
#define ll long long 

const int N=205,H=1e9+7;
int fac[N],ifac[N],inv[N],k,S[N][N],C[N][N];
ll l,r;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

struct Nod{
  int a,b;
  Nod operator +(const Nod &x) const{return {adc(a,x.a),adc(b,x.b)};}
  Nod operator +(const int &x) const{return {adc(a,x),b};}
  Nod operator -(const Nod &x) const{return {dec(a,x.a),dec(b,x.b)};}
  Nod operator -(const int &x) const{return {dec(a,x),b};}
  Nod operator *(const Nod &x) const{return {adc(mul(a,x.a),mul(5,mul(b,x.b))),adc(mul(a,x.b),mul(b,x.a))};}
  Nod operator *(const int &x) const{return {mul(a,x),mul(b,x)};}
  Nod operator /(const Nod &x) const{
	int res=qpow(dec(mul(x.a,x.a),mul(5,mul(x.b,x.b))));
    return {mul(res,dec(mul(a,x.a),mul(5,mul(b,x.b)))),mul(res,dec(mul(b,x.a),mul(a,x.b)))};
  }
};

Nod qpow(Nod a,ll b=H-2){
  Nod res={1,0};
  while(b){if(b&1) res=res*a;a=a*a,b>>=1;}
  return res;
}

Nod Sum(Nod a,ll n){return (qpow(a,n+1)-1)/(a-1);}

Nod calc(Nod x,ll l,ll r){
  if(x.a==1&&x.b==0) return x*((r-l+1)%H);
  return Sum(x,r)-Sum(x,l-1);
}

void init(){
  fac[0]=inv[1]=1;
  for(int i=1;i<N;i++) fac[i]=mul(fac[i-1],i);
  for(int i=2;i<N;i++) inv[i]=mul(inv[H%i],dec(0,H/i));
  ifac[N-1]=qpow(fac[N-1]);
  for(int i=N-1;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
  
  C[0][0]=1;
  for(int i=1;i<N;i++){
	C[i][0]=1;
	for(int j=1;j<=i;j++) C[i][j]=adc(C[i-1][j],C[i-1][j-1]);
  }
  S[0][0]=1;
  for(int i=1;i<N;i++)
    for(int j=1;j<=i;j++) S[i][j]=dec(S[i-1][j-1],mul(i-1,S[i-1][j]));
}

void SOLVE(){
  Nod a={0,inv[5]},b={0,dec(0,inv[5])},x={inv[2],inv[2]},y={inv[2],dec(0,inv[2])};
  int ans=0;
  for(int i=0;i<=k;i++){
    Nod cur={0,0};
	for(int j=0;j<=i;j++){
      Nod nw=qpow(a,j)*qpow(b,i-j)*calc(qpow(x,j)*qpow(y,i-j),l,r)*C[i][j];
	  cur=cur+nw;
	}
	cur=cur*S[k][i];
	add(ans,cur.a);
  }
  ans=mul(ans,ifac[k]);
  cout<<ans<<'\n';
}

signed main(){
  /*2024.8.1 H_W_Y CF717A Festival Organization Fibonacci + Stirling*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  init();
  cin>>k>>l>>r,l+=2,r+=2;
  SOLVE();
  return 0;
}

加强版：P5320 [BJOI2019] 勘破神机。

也可以说是双倍经验吧。

容易发现其实 \(m=2\) 的情况就是这道题，直接贺。

而关于 \(m=3\)，\(2 \nmid n\) 时显然是不合法的，于是我们设 \(f_i\) 表示长度为 \(2i\) 时的答案，有递推方程

\[f_i =2 \sum_{j=0}^{i-2} f_j + 3 f_{i-1} \]

不要忘记乘二！！！

于是就可以用类似于解 Fibonacci 的方法，用生成函数刻画得到

\[f_n = \frac {3+ \sqrt 3} 6 (2 + \sqrt3)^n + \frac {3 - \sqrt 3} 6 (2 - \sqrt 3)^n \]

就和上面一道题一样了。代码。

CF961G Partitions - Stirling 数通项公式

给出 \(n\) 个物品，每个物品有一个权值 \(w_i\)。

定义一个集合 \(S\) 的权值 \(W(S)=|S|\sum\limits_{x\in S}w_x\)。

定义一个划分的权值为 \(W'(R)=\sum\limits_{S\in R}W(S)\)。

求将 \(n\) 个物品划分成 \(k\) 个集合的所有方案的权值和。答案对 \(10^9+7\) 取模。

\(n,k\le2\times10^5\)，\(w_i\le10^9\)。

首先，容易想到去考虑每一种情况的贡献次数，也就是对于 \(w_i\)，它贡献的次数是

\[\sum_{j=1}^n j \binom {n-1} {j-1} {n-j \brace k-1} \]

怎么化简？

我们考虑把 第二类斯特林数 展开，也就是用通项公式

\[{n \brace k } = \sum_{t=0}^k \frac {(-1)^t} {t!} \frac {(k-t)^n} {(k-t)!} \]

得到

\[\begin{aligned} & \sum_{j=1}^n j \binom {n-1} {j-1} {n-j \brace k-1} \\ = & \sum_{j=1}^n j \binom {n-1} {j-1} \sum_{t=0}^{k-1} \frac {(-1)^t} {t!} \frac {(k-1-t)^{n-j}} {(k-1-t)!} \\ = & \sum_{i=0}^{k-1} \frac {(-1)^i} {i!} \frac {1} {(k-1-i)!} \sum_{j=1}^n j \binom {n-1} {j-1} (k-1-i)^{n-j} \end{aligned} \]

考虑把 \(j\) 吸收公式吸进去，也就是得到

\[\begin{aligned} & \sum_{j=1}^n j \binom {n-1} {j-1} (k-1-i)^{n-j} \\ = & \sum_{j=1}^n (j-1) \binom {n-1} {j-1} (k-1-i)^{n-j} + \sum_{j=1}^n \binom {n-1} {j-1} (k-i-1)^{n-j} \\ = & (n-1) \sum_{j=1}^n \binom{n-2} {j-2} (k-1-i)^{n-j} + (k-i)^{n-1} \\ = & (n-1) (k-i)^{n-2} + (k-i)^{n-1} \\ = & (k-i)^{n-2} (k-i-1+n) \end{aligned} \]

于是答案就是

\[\left( \sum_{i=1}^n w_i \right) \left( \sum_{i=0}^{k-1} \frac {(-1)^i} {i!} \frac {(k-i)^{n-2}} {(k-i-1)!} (k-i-1+n) \right) \]

就可以直接计算了。然而，其实后面的一块式子可以被表示成

\[{n \brace k} + (n-1) {n-1 \brace k} \]

这可能就是组合意义的解释吧。

时间复杂度 \(\mathcal O(k \log k)\)，快速幂的 \(\log\)。代码。

Code

#include <bits/stdc++.h>
using namespace std;

const int N=2e5+5,H=1e9+7;
int n,k,sum,fac[N],ifac[N],ans=0;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

void init(){
  fac[0]=1;
  for(int i=1;i<N;i++) fac[i]=mul(fac[i-1],i);
  ifac[N-1]=qpow(fac[N-1]);
  for(int i=N-1;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
}

int main(){
  /*2024.8.2 H_W_Y CF961G Partitions Stirling*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>k;init();
  for(int i=1,x;i<=n;i++) cin>>x,add(sum,x);
  if(n==1){
	if(k==1) cout<<sum<<'\n';
	else cout<<"0\n";exit(0);
  }
  for(int i=0;i<k;i++){
	int tmp=mul(mul(ifac[i],ifac[k-1-i]),mul(qpow(k-i,n-2),k-i-1+n));
	if(i&1) del(ans,tmp);
	else add(ans,tmp);
  }
  ans=mul(ans,sum);
  cout<<ans;
  return 0;
}

CF960G Bandit Blues

给你三个正整数 \(n\)，\(a\)，\(b\)，定义 \(A\) 为一个排列中是前缀最大值的数的个数，定义 \(B\) 为一个排列中是后缀最大值的数的个数，求长度为 \(n\) 的排列中满足 \(A = a\) 且 \(B = b\) 的排列个数。\(n \le 10^5\)，答案对 \(998244353\) 取模。

挺烫的。

考虑那个数既有前缀最大值又有后缀最大值：只有 \(n\)。

于是 \(n\) 相当于把序列划分成了两个区间，前面满足前缀最大值个数为 \(a-1\)，后面满足后缀最大值个数为 \(b-1\)，两者互不干扰。

我们考虑从大往小插入，设前 \(i\) 大的数构成前缀最大值 \(j\) 个的方案数，那么有

\[f_{i,j} = (i-1)f_{i-1,j} + f_{i-1,j-1} \]

简直和第一类斯特林数的转移一模一样，于是其实

\[f_{i,j} = {i \brack j} \]

那么我们就得到答案

\[\sum_{i=1}^n {i-1 \brack a-1} {n-i \brack b-1} \binom {n-1} {i-1} \]

直接贺第一类斯特林数列，然后冲过。代码。

Code

poly Stirling(int n,int m){
  poly F(n+1);
  for(int i=1;i<=n;i++) F[i]=inv[i];
  F=Qpow_plus(F,m);
  for(int i=0;i<=n;i++) F[i]=mul(F[i],mul(fac[i],ifac[m]));
  return F;  
}

int binom(int n,int m){
  if(n<m||m<0) return 0;
  return mul(fac[n],mul(ifac[m],ifac[n-m]));
}

int main(){
  /*2024.8.2 H_W_Y P5409 第一类斯特林数·列*/
  cin>>n>>a>>b;init(N-1);
  if(!a||!b) cout<<0<<'\n',exit(0);
  --a,--b;
  poly A=Stirling(n,a),B=Stirling(n,b);
  for(int i=1;i<=n;i++) add(ans,mul(binom(n-1,i-1),mul(A[i-1],B[n-i])));
  cout<<ans;
  return 0;
}

其实想想上面这个式子的组合意义，我们还可以把它化简。

答案就是

\[{n-1 \brack a+b-2} \binom {a+b-2} {a-1} \]

完全胜利！

但似乎还是要分治 NTT 或者第一类斯特林数行求，还是 \(\mathcal O(n \log n)\)，常数会小一点。