[学习笔记] 母函数

前言：

其实前面刷的FFT的题中已经包含的母函数的题目…
后来发现还有一些别的类型，于是系统的学习一下…

定义：

根据定义，这个序列作为函数的系数，称G(x)就是序列的母函数。和一般意义上的函数相比，母函数的功能是计数。

现在我们考虑一个问题

有两个骰子 请问 两个骰子的点数之和为6的情况有多少种…

高中课本里学习了两种基本的计数方法：加法原理和乘法原理
前者是分类 后者则为分步

加法原理：有可能是1+5=5+1或2+4=4+2或3+3=3+3，累加起来，一共五种可能。
乘法原理：第一个骰子只能为 1-5 而当第一个骰子的点数确定的时候 第二个骰子的点数也是确定的 所以情况为 1*5=5

现在我们把问题的要求变一下

求n个筛子点数之和为m的方案数…

总不能用上面的方法慢慢枚举吧…

那么先从简单情形考虑，两个骰子掷出n点，有多少种可能。这相当于把n拆分成两个数的和，这时候对n枚举就很复杂，所以对两个骰子枚举。第一个骰子，6种可能，相互之间是“或”的关系，对应的是加法，但是不能直接相加，因为这无法反映两个骰子的叠加过程。

我们需要逐个分析一个骰子的不同情况，对于一个骰子，假如掷出2个点，就可以看作一个分步策略，相当于——先掷一个点，再掷一个点，所以是(●)^2，当然这样不习惯，于是用x^2来表示。那四个点是x^4,因此可以用指数对应点数。

这样就可以用（x+x^2+…+x^6）表示一个骰子的投掷过程，对于第二个骰子也是这样，最后两式相乘，x^6前面的系数就是有多少种出现6点的方案。

设 $G(x)=x^1+x^2+x^3+x^4+x^5+x^6$

$H(x)=G(x)^2$

那么最后 H(x) 的 6次方的系数 即为我们要求的答案

上面的那个问题答案 即为 $G(x)^n$里面m次方的系数

这可以看出来，对于这个多项式，我们并不关心它的值，现在关心的却是它的系数了。

母函数是一个计数工具，x的取值我们不关心，只是个占位置的东西.

类型：

母函数分为普通型母函数和指数型母函数

一些简单的题目：

普通型母函数：

HDU 1028

对于每个数能够组合出的数构造一个母函数
然后把这些母函数相乘即为答案
即$(x^0+x^1+x^2+x^3+...x^n)*(x^0+x^2+x^4+...)*...*(x^0+x^n)$
N比较小 N^3暴力乘即可
FFT貌似不行，因为长度问题…

#include <cstdio>
#include <iostream>
#include <cstring>
const int maxm=1e5+100; 
int A[maxm],B[maxm];
int main()
{
    int n,m;
    while((scanf("%d",&n))!=EOF)
    {
        for(int i=0;i<=n;i++) A[i]=1,B[i]=0;
        for(int i=2;i<=n;i++)
        {
            for(int j=0;j<=n;j++)
             for(int k=0;k<=n;k+=i)
              B[j+k]+=A[j];
            for(int j=0;j<=n;j++) A[j]=B[j],B[j]=0;
        }
        printf("%d\n",A[n]);
    }
}

HDU 1398

与上题基本一致，只是能用的数只有完全平方数

#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
const int maxm=1e5+100; 
int A[maxm],B[maxm];
int main()
{
    int n,m;
    while((scanf("%d",&n))&&n)
    {
        for(int i=0;i<=n;i++) A[i]=1,B[i]=0;
        for(int i=2;i<=n;i++)
        {
            if(std::sqrt(i)!=std::ceil(std::sqrt(i))) continue;
            for(int j=0;j<=n;j++)
             for(int k=0;k<=n;k+=i)
              B[j+k]+=A[j];
            for(int j=0;j<=n;j++) A[j]=B[j],B[j]=0;
        }
        printf("%d\n",A[n]);
    }
}

HDU 1085

构造三个多项式
$(x^0+x^1+x^2+...x^a)*(x^0+x^{2*1}+x^{2*2}+...+x^{2*b})**(x^0+x^{2*1}+x^{2*2}+...+x^{2*c})$
然后暴力找系数为0的即可
FFT常数太大…

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>
#define double long double
const double PI=std::acos(-1);
const int maxm=3*5010*2;
struct complex{
    double real,imag;
    complex(){};
    complex(double _real,double _imag):real(_real),imag(_imag){}
}; 
inline complex operator + (complex x,complex y){return (complex){x.real+y.real,x.imag+y.imag};}
inline complex operator - (complex x,complex y){return (complex){x.real-y.real,x.imag-y.imag};}
inline complex operator * (complex x,complex y){return (complex){x.real*y.real-x.imag*y.imag,x.real*y.imag+y.real*x.imag};}
complex A[maxm],B[maxm],C[maxm];
int rev[maxm],len;
inline void FFT(complex *a,int n,int f)
{
    for(int i=0;i<n;i++) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1)
    {
        complex wn=(complex){std::cos(PI/i),f*std::sin(PI/i)};
        for(int j=0;j<n;j+=(i<<1))
        {
            complex w=(complex){1,0};
            for(int k=0;k<i;k++,w=w*wn)
            {
                complex x=a[j+k];
                complex y=a[i+j+k]*w;
                a[j+k]=x+y;
                a[i+j+k]=x-y;
            }
        }
    }
    if(f==-1) for(int i=0;i<n;i++) a[i].real=a[i].real/n;
}
int numa,numb,numc;
int m;
inline void work()
{
    memset(A,0,sizeof(A));
    memset(B,0,sizeof(B));
    memset(C,0,sizeof(C));
    for(int i=0;i<=numa;i++) A[i*1].real=1;
    FFT(A,m,1);
    for(int i=0;i<=numb;i++) B[i*2].real=1;
    FFT(B,m,1);
    for(int i=0;i<=numc;i++) C[i*5].real=1;
    FFT(C,m,1);
    for(int i=0;i<m;i++) B[i]=B[i]*C[i],A[i]=A[i]*B[i];
    FFT(A,m,-1);
    for(int i=1;i<m;i++) 
    if((int)round(A[i].real)==0)
    {
        printf("%d\n",i);
        return;
    } 
}
int main()
{

    for(m=1;m<=3*5001;m<<=1) len++;
    for(int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
    while((scanf("%d%d%d",&numa,&numb,&numc))&&(numa||numb||numc)) work();
    return 0;
}

指数型母函数：

形式基本为
$g(x)=\sum_{0}^{inf} \frac{x^n}{n!}=a_0+a_1*x^1+a_2*\frac{x^2}{2!}*...$

POJ 3734

直接搞
$(x^0+x^1+\frac{x^2}{2!}+...)^2*(x^0+\frac{x^2}{2!}+\frac{x^4}{4!}+...)^2=e^{2x}+(\frac{e^x+e^{-x}}{2})^2=\frac{e^{4x}+2e^{2x}+1}{4}$
$\frac{e^{4x}+2e^{2x}+1}{4}=\frac{1}{4}+\sum_{0}^{inf}\frac{4^n+2^{n+1}}{4}\frac{x^n}{n}$
最后的答案即为第n项的系数 * n!

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
const int mod=10007;
inline int fastpow(int x,int y)
{
    int ans=1;
    for(;y;y>>=1,x=(1ll*x*x)%mod) if(y&1) ans=(1ll*ans*x)%mod;
    return ans;
} 
int main()
{
    int t,n,inv=fastpow(4,mod-2);
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        printf("%d\n",(((fastpow(4,n)+fastpow(2,n+1))%mod)*inv)%mod);   
    }   
    return 0; 
}

HDU 1521

柿子就不写了，n^3暴力做就好了…

#include <cstring>
#include <algorithm>
#include <iostream>
#include <cstdio>
const int maxm=1e5+100; 
double a[maxm],b[maxm],p[maxm];
int num[maxm];
int main()
{
    p[0]=1;
    for(int i=1;i<=10;i++) p[i]=p[i-1]*(double)i;
    int n,m;
    while((scanf("%d%d",&n,&m))!=EOF)
    {
        for(int i=1;i<=n;i++) scanf("%d",&num[i]),num[i]=std::min(num[i],m);
        memset(a,0,sizeof(a)),memset(b,0,sizeof(b));
        for(int i=0;i<=num[1];i++) a[i]=1.0/p[i];
        for(int i=2;i<=n;i++)
        {
            for(int j=0;j<=m;j++)
             for(int k=0;k<=num[i];k++)
              b[j+k]+=a[j]*1.0/p[k];
            for(int j=0;j<=m;j++) a[j]=b[j],b[j]=0;
        }
        printf("%.0lf\n",a[m]*p[m]);
    }
    return 0;
}