【x义x讲坛】生成函数入门

· · 个人记录

\color{green}\Large\texttt{『菜鸡的blog』}

声明:该博客使用的例题大多来自《信息学奥赛数学一本通》。

众所周知:

(x+1)^N=C_N^0+C_N^1x+C_N^2x^2+\dots+C_N^{N-1}x^{N-1}+C_N^Nx^N

然后,你可能在高中作业里做到让你求C_N^0-C_N^1+C_N^2-C_N^3+\dots的题目。又是非常众所周知地,你把x=-1带进上面那个柿子你就会发现,等式左边是0,右边是你要求的柿子,于是就瞬间得出答案了。

扩展一下,我们可以用一个函数来代表一个数列,然后把对数列的乱搞变成对这个函数的代数运算,比如加法,减法,乘法,求导,积分之类的。于是可以随便选一批“标志函数”f_k(x),然后我们定义:

F(x)=a_0+a_1f_1(x)+a_2f_2(x)+\dots

是数列\{a_0,a_1,\dots\}(可以是无穷的)的母函数。比如,当标志函数f_k(x)=x^k的时候,(x+1)^N就是数列\{C_N^0,C_N^1,\dots,C_N^N\}的母函数。

§ 1 § 普通型母函数

普通型母函数的特色是标志函数f_k(x)=x^k。比如上面的那个例子就是。

§ 1.1 § 一些有特色的数列和它们的普通型母函数

考虑一个非常sb的数列\{1,1,1,\dots\}。它的普通型母函数是:

F(x)=1+x+x^2+x^3+\dots

显然,

xF(x)=x+x^2+x^3+\dots

两个减一下得到

(1-x)F(x)=1,F(x)=\dfrac{1}{1-x}

这个结果在某种意义上没有问题,因为如果你把一个-1<x<1x带进去的话确实正确。但是如果你把,比如x=2带进去的话会得出非常诡异的结果:1+2+4+8+\dots=-1。事实上,你不应该把x=2带进去,因为原来那个东西在x=2的时候根本不收敛。

但是,一般来讲,在做题的时候你是不用考虑这些的,因为我们主要想要的是x^k前面的那个系数而不是函数值本身。

再列出一些:

\boxed{\begin{cases}\dfrac{1}{1-x}\quad\quad=1+x+x^2+x^3+\dots\\\dfrac{1}{(1-x)^{N+1}}\ \ =\sum_{i=0}^\infty C_{N+i}^Nx^i\\\dfrac{1-x^{N+1}}{1-x}\quad\ =1+x+x^2+\dots+x^N\\\dfrac{1}{1-x^2}\quad\ \ =1+x^2+x^4+\dots\\\dfrac{x}{1-x^2}\quad\ \ =x+x^3+x^5+\dots\end{cases}}

§ 1.2 § 普通型母函数的乘法

普通型母函数的乘法和“多重集的组合问题”有关。不懂什么是“多重集的组合问题”没关系,我也不懂。我们看一道例题了解一下吧。

【例一】

现在有重量为1g的砝码两个,3g的砝码两个,5g的砝码一个。求:

【例一解答】

先放结论:

(1+x+x^2)(1+x^3+x^6)(1+x^5)

的展开式中,项数即为第一小问,x^6的系数即为第二小问。请读者自行验证。

为什么?

这可以看做是一个动态规划问题,设当前称出k克的方案数为a_k,那么如果新考虑一种砝码,质量为b,有c个,那么我可以选择不加这种砝码,加一个这种砝码,加两个……加c个,也就是把a_k贡献给新的a_k,a_{k+b},a_{k+2b},\dots,a_{k+bc}。这恰好对应了数列\{a_i\}的普通型母函数乘以1+x^b+x^{2b}+\dots+x^{bc}

其他题目也有类似的思想,此处不表。

§ 1.3 § 其他

有时候,普通型母函数只是刚好“凑到”了某个数列的递推式而已。

【例二】

求在十进制下,N位的(不允许有前导0),写出来有偶数个5的数字的个数。

【例二解答】

设答案为a_n,把“偶数”改成“奇数”的话答案为b_n。显然地a_1=8,b_1=1

对于一个n位的,有偶数个5的数,它只有两种可能:

也就是说,

a_n=9a_{n-1}+b_{n-1}

对称地,

b_n=9b_{n-1}+a_{n-1}

如果我们构造a,b的普通型母函数,那么通过A(x)xA(x)我们就能把a_na_{n-1}“对齐”。也就是说:

A(x)-9xA(x)-xB(x)=a_1=8 B(x)-9xB(x)-xA(x)=b_1=1

然后一顿解,解得

A(x)=\dfrac{-71x+8}{(1-8x)(1-10x)} B(x)=\dfrac{1-x}{(1-8x)(1-10x)}

然而这玩意解出来看似毛用没有,因为我们实际上需要的是x^n的系数而不是现在这个鬼玩意。但是我们可以把它化成:

A(x)=\dfrac 1 2(\dfrac 7 {1-8x}+\dfrac 9 {1-10x})

观察一下\dfrac 1 {1-8x}。如果把8x看成x,它显然就是我们熟悉的1+x+x^2+x^3+\dots。于是其实\dfrac 7 {1-8x}就是7*(1+8x+(8x)^2+(8x)^3+\dots)

于是,

A(x)=\dfrac 1 2\sum_{i=0}^\infty 7*(8x)^i+9*(10x)^i a_n=\dfrac 1 2(7*8^n+9*10^n)

类似地:

【例三】

求斐波那契数列:

\begin{cases}f[0]=1\\f[1]=1\\f[n]=f[n-1]+f[n-2]\quad\quad(\text{其他情况})\end{cases}

的通项公式。

【例三解答】

类似上面的方法,设F(x)=f[0]+f[1]x+f[2]x^2+\dots

F(x)-xF(x)-x^2F(x)=f[0]=1 F(x)=\dfrac 1 {1-x-x^2}

类似地,我们也要把这玩意化开。设:

F(x)=A*(\dfrac 1 {x-a}-\dfrac 1 {x-b})

则易解得

\begin{cases}a=\dfrac {1+\sqrt 5} 2\\b=\dfrac {1-\sqrt 5} 2\\A=\dfrac {\sqrt 5} 5\end{cases}

然后就很简单了吧。

你可能听过斐波那契数列通项公式的特征方程解法,反正当时我听这个方法的时候认为这简直就是玄学,这两天才终于发现这个所谓特征方程方法就是生成函数里的一个小结论……

另外,容易发现把一个数列的普通型母函数乘以1+x+x^2+x^3+\dots相当于求其前缀和。请读者自行验证。类似地,乘以1-x则相当于求其差分。所以,你应该已经可以独立完成P5488了。不就是加个NTT吗,有什么难的!

P5488代码如下:

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

const int p=1004535809,Len=262144;
const int g=3,ig=334845270;

int read(){
    int num=0,neg=1;char c=getchar();
    while(!isdigit(c)){if(c=='-') neg=-1;c=getchar();}
    while(isdigit(c)) num=num*10+c-'0',c=getchar();
    return num*neg;
}
int read_p(){
    int num=0;bool neg=0;char c=getchar();
    while(!isdigit(c)){if(c=='-') neg=1;c=getchar();}
    while(isdigit(c)) num=(1LL*num*10%p+c-'0')%p,c=getchar();
    return neg?(p-num):num;
}

int inv[300005];
int qpow(int a,int k){
    int ans=1;
    while(k){
        if(k&1) ans=1LL*ans*a%p;
        a=1LL*a*a%p;
        k>>=1;
    }
    return ans;
}

int N,K,X=1,len;bool flg;
int A[300005];
int B[300005];
int R[300005];

void NTT(int a[],bool flg){
    for(int i=0;i<X;i++)
        if(i<R[i]) swap(a[i],a[R[i]]);
    for(int mid=1;mid<X;mid<<=1){
        int Wn=qpow((flg?ig:g),(p-1)/(mid<<1));
        for(int j=0;j<X;j+=(mid<<1)){
            int w=1;
            for(int k=0;k<mid;k++){
                int a0=a[j+k],a1=1LL*w*a[j+mid+k]%p;
                a[j+k]=(a0+a1)%p;a[j+mid+k]=(a0-a1+p)%p;
                w=1LL*w*Wn%p;
            }
        }
    }
}

int main(){
    inv[1]=1;
    for(int i=2;i<=Len;i++)
        inv[i]=1LL*(p-p/i)*inv[p%i]%p;

    N=read(),K=read_p(),flg=read();
    while(X<(N<<1)) X<<=1,len++;
    for(int i=0;i<X;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(len-1));
    for(int i=0;i<N;i++)
        A[i]=read();
    B[0]=1;
    if(flg){
        for(int i=1;i<N;i++)
            B[i]=1LL*B[i-1]*(K-i+1+p)%p*inv[i]%p;
        for(int i=1;i<N;i+=2)
            B[i]=p-B[i];
    }
    else{
        for(int i=1;i<N;i++)
            B[i]=1LL*B[i-1]*(i+K-1)%p*inv[i]%p;
    }

    NTT(A,0);NTT(B,0);
    for(int i=0;i<X;i++) A[i]=1LL*A[i]*B[i]%p;
    NTT(A,1);
    int d=qpow(X,p-2);
    for(int i=0;i<N;i++)
        printf("%d ",1LL*A[i]*d%p);
}

§ 2 § 指数型母函数

指数型母函数的标志函数是f_k(x)=\dfrac {x^k} {k!}。(我们认为0的阶乘是1。)

§ 2.1 § 一些有特色的数列和它们的指数型母函数

仍然考虑sb数列\{1,1,1,\dots\}。它的指数型母函数是:

F(x)=1+x+\dfrac {x^2}{2!}+\dfrac{x^3}{3!}+\dfrac{x^4}{4!}+\dots

F(x)求导可以得到:

F'(x)=1+x+\dfrac {x^2}{2!}+\dfrac{x^3}{3!}+\dfrac{x^4}{4!}+\dots=F(x)

又由于F(0)=1,所以解得F(x)=e^x。这也是为什么这类母函数被称为“指数”型母函数的原因。

再列出一些。

\boxed{\begin{cases}e^x\quad\quad\quad\ \ =1+x+\dfrac{x^2}{2!}+\dfrac{x^3}{3!}+\dots\\\dfrac{e^x+e^{-x}}2\ \ =1+\dfrac{x^2}{2!}+\dfrac{x^4}{4!}+\dfrac{x^6}{6!}+\dots\\\dfrac{e^x-e^{-x}}2\ \ =x+\dfrac{x^3}{3!}+\dfrac{x^5}{5!}+\dfrac{x^7}{7!}+\dots\end{cases}}

§ 2.2 § 指数型母函数的乘法

指数型母函数的乘法和“多重集的排列问题”有关。

【例四】

用两个1,两个3,一个5:

【例四解答】

展开(1+x+\dfrac{x^2}{2!})(1+x+\dfrac{x^2}{2!})(1+x)\dfrac{x^4}{4!}的系数就是答案。

一点点简单的排列组合知识即可解释,于是不再赘述。

【例五】

求由ACTG组成的n位字符串的个数,其中A和C必须出现偶数次。

【例五解答】

A和C可以用之前的\dfrac{e^x+e^{-x}}2表示,于是得出生成函数为

\dfrac{(e^x+e^{-x})^2e^{2x}}4

括号拆成

\dfrac{e^{4x}+2e^{2x}+1}{4}

还是还原回sigma的形式:

\dfrac 1 4+\sum_{i=0}^\infty \dfrac{4^i+2*2^i}4*\dfrac{x^i}{i!}

于是

a_n=\begin{cases}\dfrac{4^i+2*2^i} 4\quad\quad(i\neq0)\\1\quad\quad\quad\quad\quad\quad\ \!(i=0)\end{cases}
【例六】

UVA1305 Chocolate

【例六解答】

由于每一个抽取排列都是等概率的,所以我们把抽n个剩m个的方案数除以c^n就是概率。

把问题转化一下:其实就是m个颜色抽了奇数次,c-m个颜色抽了偶数次。于是设计出母函数:

\dfrac{(e^x-e^{-x})^m*(e^x+e^{-x})^{c-m}}{2^c}

然后二项式定理展开一下,预处理一下组合数就可以做了。

不幸的是作者代码的精度死活卡不过去所以并没有(在UVA上)AC……(POJ上倒是A了)

给出不能AC的代码:

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

int C,N,M;
long long CC[105][105];
db ANS[105][105];
db ans;

double qpow(double a,int k){
    double tmp=1;
    while(k){
        if(k&1) tmp=tmp*a;
        a=a*a;
        k>>=1;
    }
    return tmp;
}

int main(){
    freopen("choc.in","r",stdin);
    freopen("choc.out","w",stdout);
    CC[0][0]=1;
    for(int i=1;i<=100;i++){
        CC[i][0]=1;
        for(int j=1;j<=i;j++)
            CC[i][j]=CC[i-1][j-1]+CC[i-1][j];
    }

    scanf("%d",&C);
    while(C){
        scanf("%d%d",&N,&M);
        if(N==0&&M==0){printf("1.000\n");scanf("%d",&C);continue;}
        if(M>C||M>N){printf("0.000\n");scanf("%d",&C);continue;}
        ans=0.0;
        for(int k=0;k<=C-M;k++)
        for(int i=0;i<=M;i++)
            if(i&1) ans-=qpow(1-2.0*(db)(k+i)/C,N)*CC[C-M][k]*CC[M][i];
            else ans+=qpow(1-2.0*(db)(k+i)/C,N)*CC[C-M][k]*CC[M][i];
        printf("%.3Lf\n",fabs(CC[C][M]/qpow(2,C)*ans));
        scanf("%d",&C);
    }

    return 0;
}

§ 3 § 小结

于是在作者能力范畴内的内容就讲完了……wtcl

由于这些母函数很多时候都是多项式,所以经常和多项式那一堆板子扯上关系,而作者连FFT都不会……

先写到这里算了。

UPD:更多内容会在生成函数再入门更新。