组合数学乱学

0.前言与ntt

组合数学过于差得我，必须重修！

ntt板子？其中MOD，G,INVG需要自行定义

namespace Poly{
    
    int rev[MN];

    int ksm(int a,int b){
        int ret=1;
        while(b){
            if(b&1) ret=ret*a%MOD;
            a=a*a%MOD;
            b>>=1;
        }
        return ret;
    }

    void dorev(int f[],int len){
        for(int i=0;i<len;i++){
            rev[i]=rev[i>>1]>>1;
            if(i&1){
                rev[i]|=len>>1;
            }
        }
        for(int i=0;i<len;i++){
            if(i<rev[i]) swap(f[i],f[rev[i]]);
        }
    }

    void NTT(int f[],int len,int op){
        dorev(f,len);
        for(int i=1;i<len;i<<=1){
            int Wn=ksm((op==1?G:INVG),(MOD-1)/(i<<1));
            for(int j=0;j<len;j+=(i<<1)){
                int w=1;
                for(int k=0;k<i;k++,w=(w*Wn)%MOD){
                    int x=f[j+k],y=w*f[j+k+i]%MOD;
                    f[j+k]=(x+y)%MOD;
                    f[j+k+i]=(x-y+MOD)%MOD;
                }
            }
        }
        if(op==-1){
            int invlen=ksm(len,MOD-2);
            for(int i=0;i<len;i++) f[i]=f[i]*invlen%MOD;
        }
    }

    // f is the output array
    void Mul(int f[],int g[],int n,int m){
        m+=n;
        n=1;
        while(n<=m) n<<=1;
        NTT(f,n,1);
        NTT(g,n,1);
        for(int i=0;i<=n;i++) f[i]=f[i]*g[i]%MOD;
        NTT(f,n,-1);
    }
}

1. 二项式反演

定义与形式

二项式反演用于解决“某个物品恰好若干个”这类计数问题。

两种形式，我们统计$f_{n}$ 为恰好$n$ 种的方案数量：

• 至多：$g_{n}=\sum\limits_{i=0} ^n \binom{n}{i} f_{i} \Leftrightarrow f_{n}=\sum\limits_{i=0}^n (-1)^{n-i} \binom{n}{i} g_{i}$
• 至少：$g_{k}=\sum\limits_{i=k} ^n \binom{i}{k} f_{i} \Leftrightarrow f_{k}=\sum\limits_{i=k}^n (-1)^{i-k} \binom{i}{k} g_{i}$

速记：至少往上走，至多往下走，容斥系数就是组合数上下相减。

有的时候，题目不会明确标出恰好，需要自行探索。

至少和至多的选取根据题目类型来选取，哪个简单那个就可以。

例题

HDU1465

求长为$n$ 的序列的错排的数量，其中$1\le n \le 20$.

设$f_n$ 表示恰好有$n$ 个位置错开，$g_{n}$ 表示至多$n$ 个位置错开，注意到，$g_n$ 可以瞎排列即可，所以$g_{n}=n!$。

答案即为：

$$f_{n}=\sum\limits_{i=0}^n \binom{n}{i} (-1)^{n-i} i!$$

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=25;
int ans,pw[MN],n;

int getC(int a,int b){
    if(a<b) return 0;
    return pw[a]/(pw[b]*pw[a-b]);
}

void init(){
    pw[0]=1;
    for(int i=1;i<MN;i++){
        pw[i]=pw[i-1]*i;
    }
}

void solve(){
    ans=0;
    for(int i=0;i<=n;i++){
        ans+=getC(n,i)*((n-i)&1?-1:1)*pw[i];
    }
    cout<<ans<<'\n';
}

signed main(){
    init();
    while(cin>>n){
        solve();
    }
    return 0;
}

洛谷P4859 已经没有什么好害怕的了

设$a_{i}>b_{i}$ 的数量为$x$ 个，所以有$x+(x-k)=n$

即$x=\dfrac{n+k}{2}$，也就是说$a_{i}>b_{i}$ 的对数恰好为$k$ 个。

先两个序列排序，双指针预处理$a_{i}> b_i$ 的个数$c_{i}$

考虑动态规划，我们设$g(i,j)$ 表示前$i$ 个数，强制钦定$j$ 个数使得$a_{k}>b_{k}$，而其他位置任意，即至少满足形式。

分类讨论得到转移方程：

$$g(i,j)=g(i-1,j)+g(i-1,j-1)\times(c_{i}-(j-1))$$

那么：

$$g_{i}=g(n,i) \times (n-i)!$$

至少与恰好，二项式反演即可，故代码：

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=4090,MOD=1e9+9;
int n,K,ans,fc[MN],f[MN][MN],a[MN],b[MN],c[MN];

int ksm(int a,int b){
    int ret=1;
    while(b){
        if(b&1) ret=ret*a%MOD;
        a=a*a%MOD;
        b>>=1;
    }
    return ret;
}

void init(){
    fc[0]=1;
    for(int i=1;i<MN;i++){
        fc[i]=fc[i-1]*i%MOD;
    }
}

int getC(int a,int b){
    return fc[a]*ksm(fc[b],MOD-2)%MOD*ksm(fc[a-b],MOD-2)%MOD;
}

signed main(){
    init();
    cin>>n>>K;
    K=(n+K)/2;
    for(int i=1;i<=n;i++){
        cin>>a[i];
    }
    for(int i=1;i<=n;i++){
        cin>>b[i];
    }
    sort(a+1,a+1+n);
    sort(b+1,b+1+n);
    for(int i=1,p=1;i<=n;i++){
        while(p<=n&&b[p]<a[i]){
            p++;
        }
        c[i]=(p-1);
    }
    for(int i=0;i<=n;i++) f[i][0]=1;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=i;j++){
            f[i][j]=(f[i][j]+f[i-1][j-1]*max(c[i]-j+1,0ll)%MOD+f[i-1][j])%MOD;
        }
    }
    for(int i=K;i<=n;i++){
        int tmp=f[n][i]*fc[n-i]%MOD; 
        ans=(ans+tmp*getC(i,K)*((i-K)&1?-1:1)%MOD)%MOD;
    }
    cout<<(ans+MOD)%MOD;
    return 0;
}

BZOJ2839 组合计数

一个有$N$ 个元素的集合有$2^N$ 个不同子集（包含空集），现在要在这$2^N$ 个集合中取出若干集合（至少一个），使得它们的交集的元素个数为$K$，求取法的方案数，答案模$10^9+7$。

设交集元素个数恰好为$k$ 个元素的方案数为$f_{k}$，观察至多不好求，考虑设至少$g_{k}$

对于$g_{k}$ 只需要强制钦定谁为交集元素，首先$\binom{n}{k}$ 飞上去，让后其他是任意选择的，数量为$2^{n-k}$。我们要求的是方案数，对于这是一个族，因为集合选择是相互独立的，但是至少要选取一个集合不能为空集，所以方案数即为$2^{2^{n-k}}-1$ 个集族，乘法原理：

$$g_{k}=\binom{n}{k} (2^{2^{n-k}}-1)$$

二项式反演即可：

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=1e6+15,MOD=1e9+7;
int n,K,pw[MN],inv[MN],g[MN];

int ksm(int a,int b,int P){
    int ret=1;
    while(b){
        if(b&1) ret=ret*a%P;
        a=a*a%P;
        b>>=1;
    }
    return ret;
}

void init(int n){
    pw[0]=1;
    for(int i=1;i<=n;i++) pw[i]=pw[i-1]*i%MOD;
    inv[n]=ksm(pw[n],MOD-2,MOD);
    for(int i=n;i>=1;i--){
        inv[i-1]=inv[i]*i%MOD;
    }
}

int getC(int a,int b){
    if(a<b) return 0;
    return pw[a]*inv[b]%MOD*inv[a-b]%MOD;
}

namespace ly
{
    namespace IO
    {
        #ifndef LOCAL
            constexpr auto maxn=1<<20;
            char in[maxn],out[maxn],*p1=in,*p2=in,*p3=out;
            #define getchar() (p1==p2&&(p2=(p1=in)+fread(in,1,maxn,stdin),p1==p2)?EOF:*p1++)
            #define flush() (fwrite(out,1,p3-out,stdout))
            #define putchar(x) (p3==out+maxn&&(flush(),p3=out),*p3++=(x))
            class Flush{public:~Flush(){flush();}}_;
        #endif
        namespace usr
        {
            template<typename type>
            inline type read(type &x)
            {
                x=0;bool flag(0);char ch=getchar();
                while(!isdigit(ch)) flag^=ch=='-',ch=getchar();
                while(isdigit(ch)) x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
                return flag?x=-x:x;
            }
            template<typename type>
            inline void write(type x)
            {
                x<0?x=-x,putchar('-'):0;
                static short Stack[50],top(0);
                do Stack[++top]=x%10,x/=10;while(x);
                while(top) putchar(Stack[top--]|48);
            }
            inline char read(char &x){do x=getchar();while(isspace(x));return x;}
            inline char write(const char &x){return putchar(x);}
            inline void read(char *x){static char ch;read(ch);do *(x++)=ch;while(!isspace(ch=getchar())&&~ch);}
            template<typename type>inline void write(type *x){while(*x)putchar(*(x++));}
            inline void read(string &x){static char ch;read(ch),x.clear();do x+=ch;while(!isspace(ch=getchar())&&~ch);}
            inline void write(const string &x){for(int i=0,len=x.length();i<len;++i)putchar(x[i]);}
            template<typename type,typename...T>inline void read(type &x,T&...y){read(x),read(y...);}
            template<typename type,typename...T>
            inline void write(const type &x,const T&...y){write(x),putchar(' '),write(y...),sizeof...(y)^1?0:putchar('\n');}
            template<typename type>
            inline void put(const type &x,bool flag=1){write(x),flag?putchar('\n'):putchar(' ');}
        }
        #ifndef LOCAL
            #undef getchar
            #undef flush
            #undef putchar
        #endif
    }using namespace IO::usr;
}using namespace ly::IO::usr;

signed main(){
    read(n,K);
    init(n);
    int ans=0;
    for(int i=K;i<=n;i++){
        int ksmm=ksm(2,ksm(2,n-i,MOD-1),MOD);
        int tmp=getC(n,i)*(ksmm-1)%MOD;
        ans=(ans+((i-K)&1?-1:1)*tmp%MOD*getC(i,K)%MOD+MOD)%MOD;
    }
    put((ans+MOD)%MOD);
    return 0;
}

HAOI2018 染色

这里至少好求，还是上面的恰好与至少。

我们考虑求解$g_{k}$，首先飞上去一个强制钦定$\binom{m}{k}$。那么还剩下$n-k\times S$ 个位置没有染色，这些位置可以看作同一类，方案是$(n-k \times S)!$。注意到我们还有排列，但是观察染色实际上是一个多重集排列的形式，那么公式即为$\dfrac{n!}{(S!)^k (n-k\times S)!}$。这样就满足了约束，但是还有$n-k \times S$ 空位置啊，也不能剩下，考虑染色方案，那么就是$(m-k)^{n-k \times S}$。

乘法原理，答案即为：

$$g_{k}=\binom{m}{k}\dfrac{n!}{(S!)^k (n-k\times S)!}(m-k)^{n-k \times S}$$

设$R=\min(m,\lfloor \dfrac{n}{S} \rfloor)$ 为真正能使用的颜色种类，代入形式二有：

$$\begin{aligned} f_{k}&=\sum\limits_{i=k}^R (-1)^{i-k} \binom{i}{k} g_{i} \\ & = \sum\limits_{i=k}^R (-1){i-k} \frac{i!}{k! (i-k)!} \\ &\text{考虑将 } k! \text{ 移项有：} \\ f_{k} \times k! & = \sum\limits_{i=k}^R (-1)^{i-k}\frac{i!}{(i-k)!} g_{i} \end{aligned}$$

设$A_x=\sum\limits_{i=0}^R g_{i}i! \cdot x^{i},B_{x}=\sum\limits_{i=0}^R = \dfrac{(-1)^i}{i!} \cdot x^i$，差卷积 NTT既得：

#include<bits/stdc++.h>
#include <cmath>
#define int long long
using namespace std;
constexpr int MN=1e7+15,MOD=1004535809,MODG=3,MODGinv=334845270;
int n,m,s,ans,f[MN],g[MN],pw[MN],rev[MN],inv[MN],w[MN];

int ksm(int a,int b){
    int ret=1;
    while(b){
        if(b&1) ret=ret*a%MOD;
        a=a*a%MOD;
        b>>=1;
    }
    return ret;
}

int getC(int a,int b){
    if(a<b) return 0;
    return pw[a]*inv[b]%MOD*inv[a-b]%MOD;
}

void init(){
    pw[0]=1;
    for(int i=1;i<MN;i++) pw[i]=pw[i-1]*i%MOD;
    inv[MN-1]=ksm(pw[MN-1],MOD-2);
    for(int i=MN-2;i>=0;i--){
        inv[i]=inv[i+1]*(i+1)%MOD;
    }
}

void dorev(int f[],int len){
    for(int i=0;i<len;i++){
        rev[i]=rev[i>>1]>>1;
        if(i&1){
            rev[i]|=len>>1;
        }
    }
    for(int i=0;i<len;i++){
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    }
}

void NTT(int f[],int len,int op){
    dorev(f,len);
    for(int i=1;i<len;i<<=1){
        int Wn=ksm((op==1?MODG:MODGinv),(MOD-1)/(i<<1));
        for(int j=0;j<len;j+=(i<<1)){
            int w=1;
            for(int k=0;k<i;k++,w=(w*Wn)%MOD){
                int x=f[j+k],y=w*f[j+k+i]%MOD;
                f[j+k]=(x+y)%MOD;
                f[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
if (op == -1) {
    int inv_len = ksm(len, MOD-2);
    for (int i = 0; i < len; i++) 
        f[i] = f[i] * inv_len % MOD;
}
}

int getG(int x){
    return getC(m,x)*pw[n]%MOD*ksm(inv[s],x)%MOD*inv[n-s*x]%MOD*ksm(m-x,n-s*x)%MOD;
}

signed main(){
    init();
    cin>>n>>m>>s;
    int lim=min(m,n/s);
    for(int i=0;i<=lim;i++){
        f[i]=getG(i)*pw[i]%MOD;
        g[i]=(i&1)?MOD-inv[i]:inv[i];
    }
    reverse(f,f+1+lim);
    int len;
    for(len=1;len<lim*2+2;len<<=1);
    NTT(f,len,1);
    NTT(g,len,1);
    for(int i=0;i<len;i++){
        f[i]=f[i]*g[i]%MOD;
    }
    NTT(f,len,-1);
    reverse(f,f+lim+1);
    for(int i=0;i<=lim;i++){
        int w;
        cin>>w;
        ans=(ans+f[i]*inv[i]%MOD*w)%MOD;
    }
    cout<<ans%MOD;
    return 0;
}

洛谷P9850
[ICPC 2021 Nanjing R] Ancient Magic Circle in Teyvat

给定$n$ 个点的完全图，其中$m$ 条边为红色边，其余边为蓝色边。
定义以下：
$f_{\text{红色}}$ 为四元组$(i,j,k,l)$，其中任意两点都有红色边连接的个数。
$f_{\text{蓝色}}$ 为四元组$(i,j,k,l)$，其中任意两点都有蓝色边连接的个数。
求$|f_{\text{红色}}-f_{\text{蓝色}}|$。
其中$1\le n \le 10^5,1\le m \le 2\times 10^5$。

赛时没想出正解（废话都没学二项式反演能做？）

发现蓝色很难受，显然可以考虑以下容斥，但是怎么容斥呢？
红色边的信息我们是有的，我们可以通过红色边来容斥。

但是这咋求啊？暴力枚举直接$O(n^4)$ 了www。

对于一张存在$j$ 条红色边的图，假设存在$g_i$ 个$i$ 条边的红色子图，而且只需要满足边颜色都是红色就可以了，那么有$\binom{j}{i}$ 种选择方法，那么我们不妨设$f_i$ 表示四元组存在$i$ 条边的红色子图个数，有下列式子：

$$g_i=\sum\limits_{j=i}\binom{j}{i} f_j$$

长的就很二项式反演：

$$f_i=\sum\limits_{j=i}\binom{j}{i} (-1)^{j-i} g_j$$

那么$|f_6-f_0|=|g_0-g_1+g_2-g_3+g_4-g_5|$。不难发现可以一个一个讨论（废话那你怎么求）

• $g_0$：不选红色边，瞎选4个点：$\dbinom{n}{4}$。
• $g_1$:选一个的方案数$\dbinom{m}{1}$，让后在确定2个端点瞎选：$\dbinom{m}{1}\dbinom{n-2}{2}$。
• $g_2$：分类讨论
  • 如果是一点连两条边，枚举公共点，让后再枚举以该端点出发的两个点，让后再瞎选一个：$(n-3)\sum\limits_{i=1}^{n}\dbinom{deg_i}{2}$，其中$deg_i$ 表示节点$i$ 的度数
  • 如果没有公共点，正南则反，就是原图任意选2个边减去有公共点的，即：$\dbinom{m}{2}-\sum\limits_{i=1}^n \dbinom{deg_i}{2}$。
• $g_3$：继续
  • 如果三元环，那就枚举剩下一个点为$(n-3)C_3$.
  • 如果是共用一个顶点，那么很简单直接枚举即可，结果$\sum\limits_{i=1}^n \binom{n}{3}$。
  • 如果是链，注意一下要把三元环的三个情况舍去，结果就是$\sum\limits_{(u,v)\in E}(deg_u-1)(deg_v-1)-3C_3$。
• $g_4$:
  • 如果四元环，那不用枚举直接$C_4$.
  • 如果三元环出来一个那就是$\sum\limits_{i=1}^n T_i(deg_i-2)$，其中$T_i$ 为$i$ 号点不同三元环的个数
• $g_5$ ：只能是两个三元环共用一条边枚举公共边即可，其中$f_5=\sum\limits_{i=\in \mathbb{C}_3}\binom{t_i}{2}$ ，其中$\mathbb{C}_3$ 表示求解三元环完成定向的边集，$t_i$ 表示覆盖带边$i$ 的不同三元环个数。

做完了，直接公式计算即可，注意瓶颈在三元环和四员化计算，不要超过$O(n^2)$：

#include <bits/stdc++.h>
#define pir pair<int, int>
#define ll long long
using namespace std;
constexpr int MN = 1e5 + 15, MM = 2e5 + 15;
struct Edge {
    int u, v;
} e[MM];
ll f0, f1, f2, f3, f4, f5;
int dg[MN],n,m,top,s[MN],id[MN];
ll cp[MN], ce[MM];
vector<int> adj[MN];
vector<pir> G[MN];

ll countthree() {
    ll ret = 0;
    for (int i = 1; i <= n; i++) {
        for (auto p : G[i]) id[p.first] = p.second;
        for (auto p : G[i]) {
            int v = p.first;
            for (auto pv : G[v]) {
                int w = pv.first;
                if (id[w]) {
                    ret++;
                    cp[i]++;
                    cp[v]++;
                    cp[w]++;
                    ce[p.second]++;
                    ce[pv.second]++;
                    ce[id[w]]++;
                }
            }
        }
        for (auto p : G[i]) id[p.first] = 0;
    }
    return ret;
}

ll countfour() {
    memset(id, 0, sizeof(id));
    ll ret = 0;
    for (int i = 1; i <= n; i++) {
        for (int v : adj[i]) {
            for (auto p : G[v]) {
                int w = p.first;
                if (dg[i] < dg[w] || (dg[i] == dg[w] && i < w)) {
                    ret += id[w];
                    if (!id[w]) s[++top] = w;
                    id[w]++;
                }
            }
        }
        for (int j = 1; j <= top; j++) id[s[j]] = 0;
        top = 0;
    }
    return ret;
}

int main() {
    cin >> n >> m;
    for (int i = 1; i <= m; i++) {
        cin >> e[i].u >> e[i].v;
        dg[e[i].u]++;
        dg[e[i].v]++;
        adj[e[i].u].push_back(e[i].v);
        adj[e[i].v].push_back(e[i].u);
    }

    for (int i = 1; i <= m; i++) {
        int u = e[i].u, v = e[i].v;
        if ((dg[u] == dg[v] && u > v) || dg[u] > dg[v])
            swap(u, v);
        G[u].push_back({v, i});
    }

    ll tri = countthree();
    ll quad = countfour();

    for (int i = 1; i <= n; i++) {
        f2+=1LL * dg[i] * (dg[i] - 1) / 2 * (n - 4);
        f3+=1LL * dg[i] * (dg[i] - 1) * (dg[i] - 2) / 6;
        f4+=1LL * cp[i] * (dg[i] - 2);
        for (auto p : G[i]) {
            int v = p.first;
            f3+=1LL * (dg[i] - 1) * (dg[v] - 1);
        }
    }

    for (int i = 1; i <= m; i++) {
        f5+=1LL * ce[i] * (ce[i] - 1) / 2;
    }

    f0=(__int128)n * (n - 1) * (n - 2) * (n - 3) / 24;
    f1=1LL * m * (n - 2) * (n - 3) / 2;
    f2+=1LL * m * (m - 1) / 2;
    f3+=tri * (n - 6);
    f4+=quad;
    cout << abs(f0 - f1 + f2 - f3 + f4 - f5);
    return 0;
}

2 斯特林数

小约定

我们有一些小概念，比较简单：

• 下降幂，即$x^{\underline{n}}=x(x-1)(x-2)\dots(x-n+1)$

例如$x^{\underline{3}}=x(x-1)(x-2)$。

• 上升幂，即$x^{\overline{n}}=x(x+1)(x+2)\dots(x+n-1)$。

• 普通幂，就是$x^n=x_1 \times x_2 \dots \times x_n$，没啥好说的。

第一类

定义与递推

我们定义，$\begin{bmatrix} n \\ m \end{bmatrix}$ 为把$n$ 个元素分成$m$ 个环的方案数，这$m$ 个环不同当且仅当环不能旋转得到另一个环。

递推：

$$\begin{bmatrix} n \\ m \end{bmatrix}=\begin{bmatrix} n-1 \\ m-1 \end{bmatrix} + (n-1) \times \begin{bmatrix} n-1 \\ m \end{bmatrix}$$

从$n-1$ 个元素推过来，两个空环肯定是不符合的。

空一个环我们考虑插入元素，插入法即可。

性质与求法

性质：

与阶乘的关系：

• $n! =\sum\limits_{i=0}^n \begin{bmatrix} n \\ i \end{bmatrix}$

本质就是置换吗。

下降幂与上升幂与普通幂的转化：

下降幂转普通幂：

$$x^{\underline{n}}=\sum\limits_{i=0}^n \begin{bmatrix} n \\ i \end{bmatrix} (-1)^{n-i} x^{i}$$

上升幂转普通幂：

$$x^{\overline{n}}=\sum\limits_{i=0}^n \begin{bmatrix}n \\ i \end{bmatrix} x^i$$

注意，上升幂转普通幂用第一类斯特林数。

证明：

如何求呢？

$$\sum\limits_{i=0} ^n \begin{bmatrix} n \\ i \end{bmatrix} x^i =\prod_{i=0}^{n-1} (x+i)=\frac{(x+n-1)!}{(x-1)!}$$

其实就是类似上面的递推式子，我们利用分治 FFT 即可$O(n \log n)$ 求解第一类斯特林数，代码如下。

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=3e5+15,MOD=167772161,MODG=3,MODGinv=55924054;
int n,rev[MN],inv[MN],pw[MN],f[MN],g[MN],a[MN],b[MN];

int ksm(int a,int b){
    int ret=1;
    while(b){
        if(b&1) ret=ret*a%MOD;
        a=a*a%MOD;
        b>>=1;
    }
    return ret;
}

void dorev(int f[],int len){
    for(int i=0;i<len;i++){
        rev[i]=rev[i>>1]>>1;
        if(i&1){
            rev[i]|=len>>1;
        }
    }
    for(int i=0;i<len;i++){
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    }
}

void NTT(int f[],int len,int op){
    dorev(f,len);
    for(int i=1;i<len;i<<=1){
        int Wn=ksm((op==1?MODG:MODGinv),(MOD-1)/(i<<1));
        for(int j=0;j<len;j+=(i<<1)){
            int w=1;
            for(int k=0;k<i;k++,w=(w*Wn)%MOD){
                int x=f[j+k],y=w*f[j+k+i]%MOD;
                f[j+k]=(x+y)%MOD;
                f[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
if (op == -1) {
    int inv_len = ksm(len, MOD-2);
    for (int i = 0; i < len; i++) 
        f[i] = f[i] * inv_len % MOD;
}
}

void Mul(int f[],int g[],int n,int m){
    m+=n,n=1;
    while(n<m) n<<=1;
    NTT(f,n,1);
    NTT(g,n,1);
    for(int i=0;i<n;i++) f[i]=f[i]*g[i]%MOD;
    NTT(f,n,-1);
}

void init(){
    pw[0]=1;
    for(int i=1;i<MN;i++) pw[i]=pw[i-1]*i%MOD;
    inv[MN-1]=ksm(pw[MN-1],MOD-2);
    for(int i=MN-2;i>=0;i--){
        inv[i]=inv[i+1]*(i+1)%MOD;
    }
}

void solve(int f[],int n){
    if(n==1){
        f[1]=1;
        return;
    }
    if(n&1){
        solve(f,n-1);
        for(int i=n;i>=1;i--){
            f[i]=(f[i-1]+f[i]*(n-1)%MOD)%MOD;
        }
        f[0]=f[0]*(n-1)%MOD;
    }else{
        int m=n/2,ret=1;
        solve(f,m);
        for(int i=0;i<=m;i++){
            a[i]=f[i]*pw[i]%MOD;
            b[i]=ret*inv[i]%MOD;
            ret=ret*m%MOD;
        }
        reverse(a,a+1+m);
        Mul(a,b,m+1,m+1);
        for(int i=0;i<=m;i++){
            g[i]=inv[i]*a[m-i]%MOD;
        }
        Mul(f,g,m+1,m+1);
        int lim=1;
        while(lim<(m+1)<<1) lim<<=1;
        for(int i=m+1;i<lim;i++) a[i]=b[i]=g[i]=0;
        for(int i=n+1;i<lim;i++) f[i]=0;
    }
}

signed main(){
    init();
    cin>>n;
    solve(f,n);
    for(int i=0;i<=n;i++) cout<<f[i]<<" ";
    return 0;
}

第二类

定义

说实话第二类才是用的最多的。

我们定义：$\begin{Bmatrix} n \\ m \end{Bmatrix}$，或记作$S(n,m)$。表示将$n$ 个元素分成$m$ 个互不区分的非空子集的方案数。

递推式子：

$$S(n,k)=S(n-1,k-1)+k \times S(n-1,k)$$

其中$S(n,0)=[n=0]$。

考虑证明，我们插入一个新元素，用两种方案：

• 将新元素单独放入一个子集，有$S(n-1,k-1)$ 种方案
• 将新元素放入一个现有的非空子集，共$k\times S(n-1,k)$ 种方案
  加法原理即可。

性质

性质如下：

$$m^n=\sum\limits_{i=0}^m S(n,i)\times i! \times \binom{m}{i}$$

当然可以写成：

$$m^n=\sum\limits_{i=0}^m S(n,i) \times m^{\underline{i}}$$

反演时写作，即普通幂转下降幂：

$$m^n=\sum\limits_{i=0}^n S(n,i) \times m^{\underline{i}}$$

可以这么理解，我们将$n$ 个不同小球丢进$m$ 个不同的盒子的方案数，我们枚举有效盒子的个数，从$m$ 个盒子选$i$ 个盒子，让后将$n$ 个小球丢进$i$ 个盒子（即第二类斯特林数）

然而事实上，这就是第二类斯特林数的生成函数。

那普通幂转上升幂呢：

$$m^n=\sum\limits_{i=0}^n (-1)^{n-i} \times S(n,i) \times m^{\overline{i}}$$

第二类斯特林显然和排列组合有关系：

$$S(n,m)=\frac{1}{m!}\sum\limits_{i=0}^m (-1)^i \binom{m}{i} (m-i)^n$$

如果空箱子的情况算进去，那么答案就是$\dfrac{m^n}{m!}$

选$i$ 个空盒子，让后把小球放进其他盒子里。

但我们的答案是有区别的盒子，要乘以$\dfrac{1}{m!}$ 去情况，

求第二类斯特林数

如何求呢，注意到上面的式子：

$$\begin{aligned} S(n,m) & =\frac{1}{m!}\sum\limits_{i=0}^m (-1)^i \binom{m}{i} (m-i)^n \\ & =\frac{1}{m!}\sum\limits_{i=0}^m (-1)^i \frac{m!}{k!(m-k)!} (m-i)^n \\ & = \sum\limits_{i=0}^m (-1)^k \times \frac{(m-k)^n}{k!(m-k)!} \end{aligned}$$

至此，我们可以实现以$O(n\log n)$ 的时间复杂度求出$S(n)$这一行的第二类斯特林数。

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=3e6+15,MOD=167772161,MODG=3,MODGinv=55924054;
int n,pw[MN],inv[MN],A[MN],B[MN];

namespace Poly{
    int rev[MN];

    int ksm(int a,int b){
        int ret=1;
        while(b){
            if(b&1) ret=ret*a%MOD;
            a=a*a%MOD;
            b>>=1;
        }
        return ret;
    }

void dorev(int f[],int len){
    for(int i=0;i<len;i++){
        rev[i]=rev[i>>1]>>1;
        if(i&1){
            rev[i]|=len>>1;
        }
    }
    for(int i=0;i<len;i++){
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    }
}

void NTT(int f[],int len,int op){
    dorev(f,len);
    for(int i=1;i<len;i<<=1){
        int Wn=ksm((op==1?MODG:MODGinv),(MOD-1)/(i<<1));
        for(int j=0;j<len;j+=(i<<1)){
            int w=1;
            for(int k=0;k<i;k++,w=(w*Wn)%MOD){
                int x=f[j+k],y=w*f[j+k+i]%MOD;
                f[j+k]=(x+y)%MOD;
                f[j+k+i]=(x-y+MOD)%MOD;
            }
        }
    }
if (op == -1) {
    int inv_len = ksm(len, MOD-2);
    for (int i = 0; i < len; i++) 
        f[i] = f[i] * inv_len % MOD;
}
}

void Mul(int f[],int g[],int n,int m){
    m+=n,n=1;
    while(n<m) n<<=1;
    NTT(f,n,1);
    NTT(g,n,1);
    for(int i=0;i<n;i++) f[i]=f[i]*g[i]%MOD;
    NTT(f,n,-1);
}

}using namespace Poly;

void init(){
    pw[0]=1;
    for(int i=1;i<MN;i++){
        pw[i]=pw[i-1]*i%MOD;
    }
    inv[MN-1]=ksm(pw[MN-1],MOD-2);
    for(int i=MN-2;i>=0;i--){
        inv[i]=inv[i+1]*(i+1)%MOD;
    }
}

signed main(){
    init();
    cin>>n;
    for(int i=0;i<=n;i++){
        A[i]=(i&1?MOD-inv[i]:inv[i]);
        B[i]=ksm(i,n)*inv[i]%MOD;
    }
    Mul(A,B,n,n);
    for(int i=0;i<=n;i++) cout<<A[i]<<" ";

    return 0;
}

与自然数幂的关系

$$\sum\limits_{i=0}^n i^k=\sum\limits_{j=0}^k S(k,j) \frac{(n+1)^{j+1}}{j+1}$$

证明：

斯特林反演

斯特林反演：

$$f(n)=\sum\limits_{k=0}^n S(n,k)\times g(k) \Leftrightarrow g(n)=\sum\limits_{k=0}^n (-1)^{n-k} \begin{bmatrix} n \\ k \end{bmatrix} f(k)$$

$$f(n)=\sum\limits_{k=0}^n (-1)^{n-k} S(n,k)\times g(k) \Leftrightarrow g(n)=\sum\limits_{k=0}^n \begin{bmatrix} n \\ k \end{bmatrix} f(k)$$

$$f(n)=\sum\limits_{i=n} S(i,n) g(i) \Leftrightarrow g(i)= \sum\limits_{i=n} (-1)^{i-n} \begin{bmatrix} i \\ n \end{bmatrix} f(i)$$

$$f(n)=\sum\limits_{i=n} (-1)^{i-n} S(i,n) g(i) \leftrightarrow g(n)=\sum\limits_{i=n} \begin{bmatrix} i \\ n \end{bmatrix} f(i)$$

例题以及应用

题型分类：

• 函数与斯特林数公式相同
  这类问题通常需要自设函数，通过发现与斯特林数的关系利用其性质求解
• 根据题意运用斯特林函数及公式
  这类问题通常隐晦地交代了需要运用斯特林函数求解，存在有效解与无效解的原式
  往往需要经验才能快速判断选择并化简原式
• 直接推式
  这类问题会直接给出包含或间接包含斯特林函数的公式，要求简化公式以得到优秀的时间复杂度
• 斯特林反演的运用
  容斥类问题，通常需要自设函数并找到与斯特林数的关系，从而化简求解过程

3. 贝尔数

贝尔数表示将基数为$n$ 的集合划分方法的数目，划分即将$S$ 划分成两两不相交的非空子集的族。例如$B_3=5$。

特别的，$B_{0}$ 为 1 因为空集有一种划分方法。

递推公式如下：

$$B_{n+1}=\sum\limits_{k=0}^n \binom{n}{k} B_{k}$$

同时因为每个贝尔数都为相邻第二类斯特林数的和，因为第二类斯特林树是把基数为$n$ 的集合乎分成$k$ 个非空集的方法数目。

即：

$$B_{n}=\sum\limits_{k=0}^n \begin{Bmatrix} n \\ k \end{Bmatrix}$$

当贝尔数到达第$15$ 项时，已经变得非常大。可以细心观察以下数据范围。

4. Lucas 定理

《初等数论》一书对卢卡斯$\texttt{(Lucas)}$定理是这么定义的$:$

设$p$为素数$,a,b\in N^* ,$并且

$$a=\sum\limits_{i=0}^k a_i p^i,b=\sum\limits_{i=0}^k b_i p^i$$

这里$0\leqslant a_i,b_i\leqslant p-1$都是整数,$i=0,1,2,…,k.$则$:$

$$C_a^b \equiv \prod\limits_{i=0}^k C_{a_i}^{b_i} \pmod{p}$$

或者写成这样$:$

$$C_a^b \equiv C_{a_k}^{b_k} \cdot C_{a_{k-1}}^{b_{k-1}} \cdot … \cdot C_{a_0}^{b_0}\pmod{p}$$

$Ps.$如果$b_i>a_i,$那么$C_{a_i}^{b_i} = 0 .$

本质上卢卡斯定理就是把组合数的计算拆到模数$p$ 的进制表示里去做。

性质

该书又介绍了$\texttt{Lucas}$定理的两个性质$:$

• 当且仅当存在$i \in \{ 0, 1, 2, …,k\}$使得$b_i > a_i$时$,C_a^b \equiv 0 \pmod{p}.$
• $C_a^b$ 为奇数的充要条件为二进制表示下$a$的每一位上的数都不小于$b$相应位上的数，即$b$ 在二进制表示下时$a$ 的子集。

应用

$\texttt{Lucas}$定理的主要用途在于在$O(\log_p a)$的时间求出$C_a^b \mod p$的结果$.$显然地$,$有一种方法$:$

对于$a,$我们预处理出$a_0,a_1,…,a_k,$对于$b$同理$.$然后我们直接算组合数并乘起来$,$注意随时取模$.$如果$p$较小且询问较多$,$可以考虑预处理组合数$.$

这听起来很容易$,$但做起来好难啊$.$
这时候我们把秦九韶请出来$,$可以将

$$\sum\limits_{i=0}^k a_i p^i$$

变成

$$(((a_kp+a_{k-1})p+a_{k-2})p+…)p+a_0$$

这样的形式$,$对于另一个同理$,$于是我们便可以层层除法$+$取模$,$将求解变成了一个递归形式$,$即$:$

$$Lucas(a,b,p)=Lucas(a/p,b/p,p) \cdot C_{a\%p}^{b\%p},$$

这里的$Lucas(a,b,p)$即要求的$C_a^b \mod p ,$

根据我们上面所提到的，组合数的计算拆成$p$ 进制下的计算，实质是这样的：

$$\begin{aligned} \dbinom{n}{m} \bmod p &= \dbinom{n/p}{m/p} \times \dbinom{n \bmod p}{m \bmod p} \bmod p \\ &= \dbinom{n/p^2}{m/p^2} \times \dbinom{n/p \bmod p}{m/p \bmod p} \times \dbinom{n \bmod p}{m \bmod p} \bmod p \\ &= \cdots \\ &= 1 \times \cdots \times \dbinom{n/p \bmod p}{m/p \bmod p} \times \dbinom{n \bmod p}{m \bmod p} \bmod p \\ \end{aligned}$$

在上述推导中，可以发现其中每一项都带上了$\pmod p$，即将每一项都限制在了$p$ 以内，即转化到了$p$ 进制。

例题

arc137d

注意到这个贡献的顺序是一个类似于前缀和的形式，如下图：

考虑如何异或只有系数为 1 的时候才能贡献到，注意到系数一定和操作次数有关，并且这个贡献顺序有点类似于格路计数，手摸不难有系数，前面的数对$a_n$ 的贡献就是：

$$a_n =\sum_{i=0} \binom{k+i-1}{i} a_{n-i}$$

其中因为是异或$\binom{k+i-1}{i}$，不难有这个是在$\pmod 2$ 意义下的，根据卢卡斯定理，当$\binom{n}{m} \bmod 2 =1$ 当且仅当$n \text{ bitand } m=m$。

用高维后缀和实现即可，时间复杂度$O(n \log n)$。

P3373 吉夫特

根据我们上面所说的，二进制模数下的组合数要想为 1 当且仅当$m$ 是$n$ 二进制表示下的子集。这道题由于值域很小，可以直接暴力枚举子集，让后就做完了，时间复杂度为$3^{\log |V|}$。

P4345 超粒子炮

原式子所求：

$$\sum_{i=0}^k \binom{n}{i} \pmod{2333}$$

注意到 2333 是质数，根据 Lucas 定理有：

$$\begin{aligned} p & = 2333 \\ \sum_{i=0}^k \binom{n}{i} & = \sum_{i=0}^k \binom{n/p}{i/p}\binom{n \bmod p}{i\bmod p} \pmod p \end{aligned}$$

注意到前面的式子类似于整除分块，枚举$i/p$ 有：

$$\binom{n/p}{0} \sum_{i=0}^{p-1} \binom{n \bmod p}{i\bmod p}+\binom{n/p}{1} \sum_{i=0}^{p-1} \binom{n \bmod p}{i\bmod p}+\dots+\binom{n/p}{k/p} \sum_{i=0}^{p-1} \binom{n \bmod p}{i\bmod p} \pmod p$$

不妨令$f(n,k)=\sum_{i=0}^k \binom{n}{i} \pmod p$，那么有：

$$f(n,k)=f(n/p,(k/p)-1) \cdot f(n\bmod p,p-1)+\binom{n/p}{k/p}\cdot f(n\bmod p,k\bmod p)$$

前面就是上面类似于整除分块的地方，后面是遗留的小块单独算。时间复杂度$O(p^2+T\log_{p}^2 n)$