生成函数学习笔记

· · 个人记录

生成函数学习笔记

终于还是逃不过要动这个东西了……自闭……

前置知识：形式幂级数、等比数列、导数、泰勒展开、微积分、广义二项式定理。

0.我为啥要学生成函数？

考纲范围内

生成函数通俗来说是一个序列的母函数，在 OI 中多用来解决组合计数有关的问题，并时常配合邪恶多项式一同食用。

然而多项式在考纲外，于是我摆烂了。

抛开多项式，生成函数推式子的部分还是可考且非常重要的，因此有学习的必要。

1. 普通生成函数 \mathrm{Ordinary\ Generating\ Function}

形如 \sum_{i=0}^{\infty}a_ix^i 的形式幂级数称为普通生成函数（OGF），常用来处理组合数问题。

例 1：有 n 种物品，第 i 种物品有 c_i 个，求本质不同的大小为 k 的物品子集个数，两个集合不同当且仅当存在一种物品在两种集合中的个数不同。

$n,k,\sum c_i\leq 10^5$：背包被橄榄了，需要想一个高明的做法。考虑多项式乘法的过程： $$(\sum\limits_{i=1}^n a_ix^i)(\sum\limits_{i=0}^mb_ix^i)=\sum\limits_{i=0}^{n+m}(\sum\limits_{j=0}^ia_jb_{i-j})x^i$$ 这个结构和我们的背包式子异曲同工，这启示我们是否可以通过设计一个多项式来替代背包的工作呢？ 答案是可以的。考虑形式幂级数 $G_i(x)=\sum_{i=1}^{c_i}x^i$ 它的每项系数对应了每个物品选择不同个数时的方案（都是一种），那么最后的答案就是 $[x^k]\prod_{i=1}^nG_i(x)$，正确性可以自己手玩。 那么这样我们就可以直接分治 NTT 做了，时间复杂度 $O(n\log n\log k)$。 --- 我们可以将普通生成函数写成封闭形式，等比数列求和即可。 常见生成函数的封闭形式，可以根据等比数列求和得到： $\sum\limits_{i=0}^{\infty}ck^ix^{ti}=\frac c{1-kx^t}

事实上，只有左边收敛的时候式子成立，但这不重要，因为对于形式幂级数我们只关心系数而不关心自变量的值到底是多少。

来一个简单的例题：化简 \sum_{i=0}^{\infty}(i+1)x^i。

设 S=1+2x+3x^2+\cdots ，

那么 xS=x+2x^2+3x^3+\cdots，

两式相减，得到 (1-x)S=1+x+x^2+x^3+\cdots=\frac1{1-x}，

因此 S=\frac1{(1-x)^2}。

再来一个稍难一点的：求斐波那契数列的生成函数的封闭形式。

还是如法炮制：设 S=1+x+2x^2+3x^3+5x^4+\cdots，

那么 xS=x+x^2+2x^3+3x^4+\cdots,x^2S=x^2+x^3+2x^4，

根据定义不难得到 S=xS+x^2S，

所以 S=\frac1{1-x-x^2}。

上面的东西看起来很炫酷，可是有什么用呢？

回到上面斐波那契数列的生成函数，假设现在我们想求它的通项公式，怎么办？

考虑到其实我们是知道形如 \frac c{1-kx} 这样的东西的通项的，能不能把他转化成这种东西？

直接因式分解，可以得到 S=\frac 1{(1-\frac{1+\sqrt5}2x)(1-\frac{1-\sqrt5}2x)}，

然后我们要把这个函数变成若干个 \frac c{1-kx} 加和的形式，设 S=\frac a{1-\frac{1+\sqrt 5}2x}+\frac b{1-\frac{1-\sqrt5}2x}，解得 a=\frac{1+\sqrt5}{2\sqrt 5},b=-\frac{1-\sqrt5}{2\sqrt5}，

代入原式，得到 S=\frac1{\sqrt5}(\frac {\frac{1+\sqrt5}2}{1-\frac{1+\sqrt 5}2x}-\frac {\frac{1-\sqrt5}2}{1-\frac{1-\sqrt5}2x})

里面的两个封闭形式都是我们熟悉的格式了，直接换成无穷级数形式得到 S=\frac1{\sqrt5}\sum_{i=0}^{\infty}[(\frac{1+\sqrt5}2)^{i+1}-(\frac{1-\sqrt5}2)^{i+1}]x^i

所以斐波那契数列的通项公式 S_n=\frac1{\sqrt5}[(\frac{1+\sqrt5}2)^{n+1}-(\frac{1-\sqrt5}2)^{n+1}]。

通过这个例子，相信大家对于生成函数封闭形式的作用会有一个直观的感受。

2.指数生成函数 \mathrm{Exponential\ Generating\ Function}

刚才的 OGF 主要用来解决组合有关的问题，EGF 则用来解决排列有关的问题。

形如 \sum_{i=0}^{\infty}a_i\frac{x^i}{i!} 的形式幂级数被称为指数生成函数（EGF）。

为什么说 EGF 可以很好的解决排列有关的问题呢？来一个例题。

例 2：有 n 种物品，第 i 种物品有 c_i 个，求本质不同的大小为 k 的有序物品子集个数，两个集合相同当且仅当每个位置选出的物品都相同。

考虑像例 1 一样，把若干个 EGF 乘起来，看看会发生什么。

显然对于第 i 种物品的 EGF G_i=\sum_{i=0}^{\infty}\frac{x^i}{i!}，那么我们发现 [x^k]\prod_{i=1}^nG_i=\sum_{\alpha_1+\alpha_2+\cdots+\alpha_n=k}\frac1{\prod_{i=1}^n\alpha_i!}。

我们回头看一眼多重集排列计数的公式 \sum_{\alpha_1+\alpha_2+\cdots+\alpha_n=k}\frac{k!}{\prod_{i=1}^n\alpha_i!}。

几乎一模一样，只差一个系数 k!。

因此我们算出 EGF 取对应系数乘上 k! 即可，这就是 EGF 和排列计数的微妙联系。

说完了 EGF 的作用，我们再来讨论一下有关 EGF 封闭形式的问题。

先从最简单的看起：\sum_{i=0}^{\infty}\frac{x^i}{i!} 的封闭形式。

把求和符号打开：1+\frac x{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots。

欸，这我熟啊，f(x)=e^x 在 x=0 处的泰勒展开。

那问题解决了，这个 EGF 的封闭形式就是 e^x。

一个难一点的：\sum_{i=0}^{\infty}\frac{x^{2i}}{(2i)!} 的封闭形式。

怎么消掉奇数次幂的项？自然联想到 (-1)^k 奇正偶负的特性，那么我们来看看 e^{-x} 是否符合我们的需要。

直接泰勒展开：e^{-x}=\sum_{i=0}^{\infty}(-1)^i\frac{x^i}{i!}。

感谢上帝，他正是我们想要的，那么我们只要二者相加就可以消掉奇数次幂的项了。

然后调整偶数次幂的系数，容易得到封闭形式为 \frac{e^x+e^{-x}}2。

再来看这样一个东西： h_n=\sum_{i=0}^n\dbinom ni f_ig_{n-i}，我们姑且叫他组合卷积。

考虑这个卷积和 EGF 的联系，拆掉这个组合数： h_n=\sum_{i=0}^n\frac{n!f_ig_{n-i}}{i!(n-i)!}=n!\sum_{i=0}^n\frac{f_i}{i!}·\frac{g_{n-i}}{(n-i)!}。

所以 \frac{h_n}{n!}=\sum_{i=0}^n\frac{f_i}{i!}·\frac{g_{n-i}}{(n-i)!}。

明显地，这是三个 EGF 之间的运算，即 H(x)=F(x)G(x)。

这启示我们，事实上很多集合划分的计数可以与 EGF 扯上关系。

看完了普通的加减乘运算，我们来看看 EGF 的 \exp 和 \ln 运算。

从 Bell 数开始，定义 Bell 数 w_n 为 n 个数的划分个数，你可以理解为它就是一行斯特林数的和。

考虑整一个 Bell 数的通项出来，我们先从递推式看起。

显然 w_0=1，对于 n>1 的情况，我们考虑枚举最后一个数所在的集合的大小，可以得到 w_n=\sum_{i=1}^n\dbinom {n-1}{i-1}w_{n-i}=\sum_{i=0}^{n-1}\dbinom {n-1}iw_{n-i-1}。

欸，右面那个东西不就是个组合卷积吗，于是我们直接整出来 Bell 数的 EGF W(x)=\sum_{i=0}^{\infty}w_i\frac{x^i}{i!}，那么原式变为 W(x)=\int e^xW(x)\mathrm dx+1，求积分是因为我们的卷积卷出来的位置是在 n-1 上的，要用积分右移一位；加一是因为 w_0 的边界。

我们来继续化简，两边同时求导：

\frac{\mathrm dW(x)}{\mathrm dx}=e^xW(x)

移项：

\frac1{W(x)}\mathrm dW(x)=e^x\mathrm dx

再对两边同时积分：

\ln W(x)=e^x+C $\ln 1=e^0+C

得到 C=-1，于是 W(x)=\exp(e^x-1)。

然后我们就可以使用多项式 \exp 来在 O(n\log n) 的时间内得到 w_0 到 w_n 了。

可惜我不会多项式 \exp，所以没有代码。

继续来看 EGF 的 \exp 和 \ln，这里是一个更加经典的例子：求 n 个点的简单无相连通图的个数。

设这个数量是 f_n，直接求显然是不好求的，我们考虑一个更简单的问题，不要求图联通怎么做。

这个我会！直接枚举每条边是否存在，设这个东西是 g_n=2^{\frac{n(n-1)}2}。

然后我们来看 f_n 和 g_n 的关系，直接算 f_n 还是不好算，我们考虑计算 f_n 的补集。

这个就好得多，至少我们可以比较轻松地找出一个递推式，考虑枚举 n 所在的联通块的大小，那么有：

把 $f_n$ 挪过去，考虑到 $\dbinom{n-1}{n-1}=g_0=1$，因此我们可以把 $f_n$ 放到卷积里面去，于是有 $g_n=\sum_{i=1}^n\dbinom{n-1}{i-1}f_ig_{n-i}

又是我们熟悉的组合卷积，直接换成 EGF，设 F(x)=\sum_{i=0}^{\infty}f_i\frac{x^i}{(i-1)!},G(x)=\sum_{i=0}^{\infty}g_{i}\frac{x^i}{i!},H(x)=\sum_{i=0}^{\infty}g_i\frac{x^i}{(i-1)!}，那么我们有 H(x)=F(x)G(x)，我们要求 F(x)，移项 F(x)=H(x)G(x)^{-1}，已经可以多项式求逆加上 NTT O(n\log n) 做了，但这看起来还不够简洁，而且我们的主角 \ln 还没有出现。

注意到事实上我们并不需要 H(x)，因为它可以被 G'(x) 替代。

我们重新设 F(x)=\sum_{i=0}^{\infty}f_i\frac{x^i}{i!}，G(x) 不变，就有 G'(x)=F'(x)G(x)。

移项得到 F'(x)=\frac1{G(x)}G'(x)。

熟悉的形式，直接两边积分，就得到了最终结果 F(x)=\ln G(x)。直接多项式 \ln 即可，复杂度 O(n\log n)。

3. 一些例题

洛谷 P2000 拯救世界

非常板板的 OGF。

注意到两个阵独立，每种石头也独立，其实就是一个容量为 n 的背包，每种物品要满足一个限制的集合计数。

那么我们直接写出十种限制的生成函数，显然不超过 k 个的限制直接等比数列求和，必须是 k 的倍数的则写成封闭形式即可。

运算过程： \frac1{1-x^6}·\frac{1-x^{10}}{1-x}·\frac{1-x^6}{1-x}·\frac1{1-x^4}·\frac{1-x^8}{1-x}·\frac1{1-x^2}·\frac{1-x^2}{1-x}·\frac1{1-x^8}·\frac1{1-x^{10}}·\frac{1-x^4}{1-x}=\frac1{(1-x)^5}

那么我们最后的答案就是 [x^n]\frac1{(1-x)^5}，直接套用广义二项式定理 \frac1{(1-x)^k}=\sum_{i=0}^{\infty}\dbinom{n+k-1}{n-1}x^i展开可以得到答案就是 \dbinom{n+4}4。

由于 n 很大，需要高精，由于 n 实在是太大，朴素的高精 O(\lg^2n) 还是过不去，需要用 NTT 优化，复杂度 O(\lg n\log\lg n)。

#include<iostream>
#include<cstdio>
#include<string>
#include<algorithm>
using namespace std;
#define int long long
const int G=3,mod=998244353;
inline int pw(int a,int b)
{
    int res=1;
    while(b)
    {
        if(b&1)
            res=res*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return res;
}
const int invG=pw(G,mod-2);
int tr[1000001<<1],len,n,m,f[1000001<<1],g[1000001<<1],ans[1000001<<1],invn;
string s;
inline void NTT(int *f,bool flag)
{

    for(register int i=0;i<n;++i)
        if(i<tr[i])
            swap(f[i],f[tr[i]]);
    for(register int p=2;p<=n;p<<=1)
    {
        int len=p>>1,tg=pw(flag? G:invG,(mod-1)/p);
        for(register int k=0;k<n;k+=p)
        {
            int buf=1;
            for(register int l=k;l<k+len;++l)
            {
                int tt=buf*f[l+len]%mod;
                f[l+len]=(f[l]+mod-tt)%mod;
                f[l]=(f[l]+tt)%mod;
                buf=buf*tg%mod;
            }
        }
    }
}
inline void mul(int c)
{
    for(int i=0;i<n;++i)
        g[i]=0;
    for(int i=0;i<len;++i)
        g[i]=s[i]^48;
    g[0]+=c;
    NTT(f,1);
    NTT(g,1);
    for(int i=0;i<n;++i)
        f[i]=f[i]*g[i]%mod;
    NTT(f,0);
    int k=0;
    for(int i=0;i<n;++i)
    {
        f[i]=f[i]*invn%mod+k;
        k=f[i]/10;
        f[i]%=10;
    }
}
signed main()
{
    cin>>s;
    reverse(s.begin(),s.end());
    len=s.length();
    f[0]=1;
    for(m=len*5,n=1;n<=m;n<<=1);
    invn=pw(n,mod-2);
    for(register int i=0;i<n;++i)
        tr[i]=(tr[i>>1]>>1)|((i&1)? n>>1:0);
    for(int i=0;i<len;++i)
        f[i]=s[i]^48;
    ++f[0];
    for(int i=2;i<=4;++i)
        mul(i);
    for(int i=n-1;i;--i)
    {
        f[i-1]+=f[i]%24*10;
        f[i]/=24;
    }
    f[0]/=24;
    int k=0;
    for(int i=0;i<n;++i)
    {
        f[i]+=k;
        k=f[i]/10;
        f[i]%=10;
    }
    while(!f[n])
        --n;
    for(int i=n;~i;--i)
        putchar(f[i]+'0');
    puts("");
    return 0;
}

CF947E Perpetual Subtraction

orz rqy。

设 p 轮之后当前数字是 i 的概率为 f_{p,i}，根据题意不难得到 f_{p+1,i}=\sum_{j=i}^n\frac{f_{p,j}}{j+1}。

这样计算当然不够快，我们来优化一下。

设 F(x)=\sum_{i=0}^{n}f_ix^i，那么有如下推导：

F_p(x)=\sum\limits_{i=0}^nf_{p,i}x^i =\sum\limits_{i=0}^nx^i\sum\limits_{j=i}^n\frac{f_{p-1,j}}{j+1} =\sum\limits_{i=0}^n\frac{f_{p-1,i}}{i+1}\sum\limits_{j=0}^ix^j =\sum\limits_{i=0}^n\frac{f_{p-1,i}}{i+1}·\frac{x^{i+1}-1}{x-1} =\frac1{x-1}\sum\limits_{i=0}^nf_{p-1,i}·\frac{x^{i+1}-1}{i+1}

观察后面的分式，我们注意到 \frac{x^{i+1}}{i+1}=\int_0^xt^i\mathrm dt,\frac1{i+1}=\int_0^1t^i\mathrm dt，结合牛顿-莱布尼茨公式可以得到：

F_p(x)=\frac1{x-1}\sum\limits_{i=0}^nf_{p-1,i}\int_1^xt^i\mathrm dt =\frac1{x-1}\int_1^xF_{p-1}(t)\mathrm dt

这个式子不好，因为我们只会做积分下界是 0 的式子，而且除以 (x-1) 看起来也很棘手，考虑换元消掉他们俩。

设 G_p(x)=F_p(x+1)，那么有：

G_p(x)=F_p(x+1)=\frac1x\int_1^{x+1}F_{p-1}(t)\mathrm dt=\frac1x\int_0^xG_{p-1}(t)\mathrm dt

啊哈！两个不好看的地方一起消失了，这样子我们直接来看 g_{p,i} 的递推式，G_{p-1}(t) 每一项直接积分可以得到 g_{p,i}x^i=\frac{g_{p-1,i}x^{i+1}}{x(i+1)}=\frac{g_{p-1,i}x^i}{i+1}，所以 g_{p,i}=\frac1{i+1}g_{p-1,i}，这显然是个线性变换，于是有 g_{p,i}=\frac1{(i+1)^{p}}g_{0,i}。

重大的进展！但遗憾的是距离正解还有很长的路要走——因为我们求出的是 g_{p,i} 的表达式，还需要在 f_{p,i} 与 g_{p,i} 之间建立起联系。

同时这个关系应该是双向的，因为我们首先要由 f_{0,i} 推出 g_{0,i}，最后还要把 g_{m,i} 还原成 f_{m,i}。

直接把我们定义的两个生成函数展开：

\sum\limits_{i=0}^ng_{p,i}x^i=\sum\limits_{i=0}^nf_{p,i}(x+1)^i=\sum\limits_{i=0}^nf_{p,i}\sum\limits_{j=0}^i\dbinom ijx^j=\sum\limits_{i=0}^n(\sum\limits_{j=i}^n\dbinom jif_{p,j})x^i

非常明了，g_{p,i}=\sum_{j=i}^n\dbinom jif_{p,j}。

但是这样没法快速算，我们把组合数拆了：

g_{p,i}=\sum\limits_{j=i}^n\frac{j!f_{p,j}}{i!(j-i)!}=\frac1{i!}\sum\limits_{j=i}^n\frac{j!f_{p,j}}{(j-i)!}

差卷积！我们的问题已经解决了一半，只剩下怎样由 g_{m,i} 推出 f_{m,i} 了。

注意到这个关系式是二项式反演的形式，我们直接反演一波：

f_{p,i}=\sum\limits_{j=i}^n(-1)^{j-i}\dbinom jig_{p,j}

还是不好看，我们再来拆了组合数：

f_{p,i}=\frac1{i!}\sum\limits_{j=i}^n\frac{(-1)^{j-i}j!g_{p,j}}{(j-i)!}

还是差卷积。至此问题解决，代码实现需要写一个 NTT，时间复杂度 O(n\log n\log m)。

#include<iostream>
#include<cstdio>
using namespace std;
#define int long long
const int G=3,mod=998244353;
inline int pw(int a,int b)
{
    int res=1;
    while(b)
    {
        if(b&1)
            res=res*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return res;
}
const int invG=pw(G,mod-2);
int n,m,p[100001],A[200001<<1],B[200001<<1],tr[200001<<1],fac[200001<<1],inv[200001<<1],g[200001<<1];
inline int read()
{
    int x=0;
    char c=getchar();
    while(c<'0'||c>'9')
        c=getchar();
    while(c>='0'&&c<='9')
    {
        x=(x<<1)+(x<<3)+(c^48);
        c=getchar();
    }
    return x;
}
inline void NTT(int *f,bool flag,int n)
{
    for(register int i=0;i<n;++i)
        if(i<tr[i])
            swap(f[i],f[tr[i]]);
    for(register int p=2;p<=n;p<<=1)
    {
        int len=p>>1,tg=pw(flag? G:invG,(mod-1)/p);
        for(register int k=0;k<n;k+=p)
        {
            int buf=1;
            for(register int l=k;l<k+len;++l)
            {
                int tt=buf*f[l+len]%mod;
                f[l+len]=(f[l]+mod-tt)%mod;
                f[l]=(f[l]+tt)%mod;
                buf=buf*tg%mod;
            }
        }
    }
    if(!flag)
    {
        int Inv=pw(n,mod-2);
        for(register int i=0;i<n;++i)
            f[i]=f[i]*Inv%mod;
    }
}
inline void solve(int *A,int *B,int n)
{
    int N;
    for(N=n<<1,n=1;n<=N;n<<=1);
    for(register int i=0;i<n;++i)
        tr[i]=(tr[i>>1]>>1)|((i&1)? n>>1:0);
    NTT(A,1,n);
    NTT(B,1,n);
    for(register int i=0;i<n;++i)
        A[i]=A[i]*B[i]%mod;
    NTT(A,0,n);
    N>>=1;
    for(register int i=N;i<n;++i)
        A[i]=B[i]=0;
}
signed main()
{
    n=read()+1,m=read()%(mod-1);
    fac[0]=inv[0]=1;
    for(register int i=0;i<n;++i)
    {
        p[i]=read();
        if(i)
        {
            fac[i]=fac[i-1]*i%mod;
            inv[i]=pw(fac[i],mod-2);
        }
    }
    for(register int i=0;i<n;++i)
    {
        A[n-i-1]=fac[i]*p[i]%mod;
        B[i]=inv[i];
    }
    solve(A,B,n);
    for(register int i=0;i<n;++i)
        g[i]=inv[i]*A[n-i-1]%mod*pw(pw(i+1,m),mod-2)%mod;
    for(register int i=0;i<n;++i)
    {
        A[n-i-1]=fac[i]*g[i]%mod;
        B[i]=i&1? mod-inv[i]:inv[i];
    }
    solve(A,B,n);
    for(register int i=0;i<n;++i)
        printf("%lld ",A[n-i-1]*inv[i]%mod);
    puts("");
    return 0;
}