zoukankan      html  css  js  c++  java
  • 【算法专题】多项式运算与生成函数

    【快速傅里叶变换】FFT

    参考:从多项式乘法到快速傅里叶变换 by miskcoo

    FFT 学习笔记 by Menci

    (一)多项式的表示法

    系数表示法:f(x)=a[n-1]*x^(n-1)+...+a[0],称为n-1次多项式。

    点值表示法:一个n-1次多项式在复数域中有n个根,即n个(x,y)可以唯一确定一个n-1次多项式。

    对于一个多项式,从其系数表示法到其点值表示法的变换称为离散傅里叶变换(DFT),反之称为傅里叶逆变换(IDFT)

    朴素的离散傅里叶变换,枚举实现的复杂度为O(n^2)。

    快速傅里叶变换是指以O(n log n)的复杂度实现IDF和IDFT的算法,常用Cooley-Tukey算法。

    (二)复数

    复数是形如a+bi的数,当b=0时为实数。

    定义一个平面为复平面,那么平面内的每个点(a,b)唯一对应一个复数a+bi,i可以理解为y轴上的单位长度,正如1是x轴上的单位长度。

    i的本质是在数轴上定义旋转变换,i是1逆时针旋转90°,那么i^2=-1。

    复数相加,遵循平行四边形定则。

    复数相乘,模长相乘,幅角相加。

    (三)单位根

    以圆点为起点,以复平面单位圆的n等分点为终点,作n个向量,设所得幅角为正且最小的向量对应的复数为ω(1,n),即n次单位根。(括号左为上标,右为下标)。

    图片来源:OI 中的 FFT by zball

    其中B点是单位根ω(1,n),逆时针依次为ω(2,n),ω(3,n)...,ω(n,n)=ω(n,0)=1。

    计算公式:ω(k,n)=cos ( 2kπ/n ) + i*sin ( 2kπ/n )

    单位根的性质:

    (1)消去:ω(2n,2k)=ω(n,k)

    (2)折半:ω(n,k+n/2)=-ω(n,k)

    将ω(n,0)~ω(n,n-1)这n个单位根作为代表n-1次多项式的n个点的横坐标,可以得到很好的性质。

    (四)快速傅里叶变换(FFT解决DFA)

    这部分因为不会操作数学公式,直接粘贴Menci博客QAQ。

    将n-1次多项式A(x)的系数奇偶分成两个多项式A1(x)和A2(x),则A(x)=A1(x^2)+x*A2(x^2)。

    对于k<n/2,有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k))

    同时,有A(ω(n,k+n/2))=A1(ω(n/2,k)) ω(n,k)*A2(ω(n/2,k))

    对于一个k次多项式,通过奇偶分项得到两个k/2次多项式,分别计算后再调用其值解决k次多项式,即分治解决。

    (五)傅里叶逆变换(IDFA)

    对于n-1次多项式,其n-1维系数向量{a0,a1...an-1}通过DFA得到点值向量{b0,b1...bn-1},反之操作称为IDFA。

    将点值向量作为系数,以单位根的倒数进行FFT,得到的每个数除以n,就是IDFA的结果。

    (六)迭代实现FFT

    对于多项式A(x),已知系数向量a[],求横坐标为ω(n,0)~ω(n,n-1)的点值向量b[]。

    将多项式奇偶分项后,对于k<n/2,有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k)),同时有A(ω(n,k+n/2))=A1(ω(n/2,k)) ω(n,k)*A2(ω(n/2,k)),分治边界是a[i]*ω(0,0)即a[i]。

    边界元素:FFT递归边界的数组排布恰好是原数组每个位置二进制反转后的数字,例如:

    原:00 01 10 11

    终:00 10 01 11

    蝴蝶操作:为了在合并时不引入新数组,进行一下操作。ω(l,k)=ω(n,n/l*k),预处理以n为底的ω[],IDFT时预处理倒数。

    t=ω(n/l*k)*a[l/2+k]

    a[l/2+k]=a[k]-t

    a[k]+=t

    (七)多项式乘法

    多项式的点值表示法易于进行乘法,因为对于fc(x)=fa(x)*fb(x),每个点x在多项式A,B中的点值相乘即可得到在多项式C中的点值。

    将n-1次多项式A和m-1次多项式B均视为n+m-2次多项式(高位补0),进行DFT后相乘再通过IDFT即可得到多项式C。

    #include<cstdio>
    #include<cstring>
    #include<cctype>
    #include<cmath>
    #include<queue>
    #include<stack>
    #include<set>
    #include<vector>
    #include<algorithm>
    #define ll long long
    #define lowbit(x) x&-x
    using namespace std;
    int read(){
        char c;int s=0,t=1;
        while(!isdigit(c=getchar()))if(c=='-')t=-1;
        do{s=s*10+c-'0';}while(isdigit(c=getchar()));
        return s*t;
    }
    int min(int a,int b){return a<b?a:b;}
    int max(int a,int b){return a<b?b:a;}
    int ab(int x){return x>0?x:-x;}
    //int MO(int x){return x>=MOD?x-MOD:x;}
    //void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;}
    /*------------------------------------------------------------*/
    const int inf=0x3f3f3f3f;
    const int maxn=300010;//2^18!!!
    const double PI=acos(-1);
    int n,m;
    struct cp{
        double x,y;
        cp(double a,double b){x=a;y=b;}
        cp(){x=y=0;};
        cp operator + (cp a){return cp(x+a.x,y+a.y);}
        cp operator - (cp a){return cp(x-a.x,y-a.y);}
        cp operator * (cp a){return cp(x*a.x-y*a.y,x*a.y+y*a.x);}
    }a[maxn],b[maxn];
    void fft(cp *a,int n,int f){
        int k=0;
        for(int i=0;i<n;i++){
            if(i>k)swap(a[i],a[k]);
            for(int j=n>>1;(k^=j)<j;j>>=1);
        }
        for(int l=2;l<=n;l<<=1){
            int m=l/2;
            cp wn(cos(2*PI*f/l),sin(2*PI*f/l));
            for(cp *p=a;p!=a+n;p+=l){
                cp w(1,0);
                for(int i=0;i<l/2;i++){
                    cp t=w*p[i+m];
                    p[i+m]=p[i]-t;
                    p[i]=p[i]+t;
                    w=w*wn;
                }
            }
        }
        if(f==-1){for(int i=0;i<n;i++)a[i].x/=n;}//
    }
    int main(){
        scanf("%d%d",&n,&m);n++;m++;
        for(int i=0;i<n;i++){int u;scanf("%d",&u);a[i]=cp(u,0);}
        for(int i=0;i<m;i++){int u;scanf("%d",&u);b[i]=cp(u,0);}
        int N=1;while(N<n+m)N*=2;
        fft(a,N,1);fft(b,N,1);
        for(int i=0;i<N;i++)a[i]=a[i]*b[i];
        fft(a,N,-1);
        for(int i=0;i<n+m-1;i++)printf("%d ",(int)(a[i].x+0.1));
        return 0;
    }
    View Code

    注意:

    1.数组空间必须是大于两个卷积数组长度和的2的幂(记为N),当n为1e5时数组空间为300000

    2.DFT后的操作一定要以N为单位,例如点值相乘。这个问题在多重fft的题目中很容易写错,多重fft必须把上一个的结果后半部分清零再继续。

    3.代码需要手写复数模板来减小常数,NTT还要预处理omega。枚举反二进制的方法是从高位开始模拟进位。

    (八)卷积

    对于函数f(n)和g(n),定义其卷积为函数(f⊗g)

    (fg)(n)=Σf(i)g(n-i),i=0~n。——形式幂级数

    卷积的形式和多项式乘法类似,(fg)的生成函数是f(n)和g(n)的生成函数的乘积。

    卷积是和为定值的形式,若差为定值,将其中一个数组反转后即可卷积。卷积中多余的部分数组为0,不影响答案。

    例题:【BZOJ】2194: 快速傅立叶之二

    【BZOJ】3527: [Zjoi2014]力 FFT

    高精度乘法:将数字从低位到高位编号0~len-1,每一位代入多项式系数,那么数字乘法就是多项式乘法,最后从低到高处理进位。

    (九)模数任意——fft合并

    将每个数字拆成a*32768+b,然后四次dft后4次idft合并。

    #include<cstdio>
    #include<cstring>
    #include<cctype>
    #include<cmath>
    #include<queue>
    #include<stack>
    #include<set>
    #include<vector>
    #include<complex>
    #include<algorithm>
    #define ll long long
    #define lowbit(x) x&-x
    using namespace std;
    int read(){
        char c;int s=0,t=1;
        while(!isdigit(c=getchar()))if(c=='-')t=-1;
        do{s=s*10+c-'0';}while(isdigit(c=getchar()));
        return s*t;
    }
    int min(int a,int b){return a<b?a:b;}
    int max(int a,int b){return a<b?b:a;}
    int ab(int x){return x>0?x:-x;}
    //int MO(int x){return x>=MOD?x-MOD:x;}
    //void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;}
    /*------------------------------------------------------------*/
    const int inf=0x3f3f3f3f,MOD=1e9+7,maxn=300010;//
    const double PI=acos(-1);
    namespace fft{
        complex<double>o[maxn],oi[maxn];
        void init(int n){
            for(int k=0;k<n;k++){o[k]=complex<double>(cos(2*PI*k/n),sin(2*PI*k/n));oi[k]=conj(o[k]);}
        }
        void transform(complex<double>*a,int n,complex<double>*o){
            int k=0;
            while((1<<k)<n)k++;
            for(int i=0;i<n;i++){
                int t=0;
                for(int j=0;j<k;j++)if(i&(1<<j))t|=(1<<(k-j-1));
                if(i<t)swap(a[i],a[t]);
            }
            for(int l=2;l<=n;l*=2){
                int m=l/2;
                for(complex<double>*p=a;p!=a+n;p+=l){
                    for(int i=0;i<m;i++){
                        complex<double>t=o[n/l*i]*p[i+m];
                        p[i+m]=p[i]-t;
                        p[i]+=t;
                    }
                }
            }
        }
        void dft(complex<double>*a,int n){transform(a,n,o);}
        void idft(complex<double>*a,int n){
            transform(a,n,oi);
            for(int i=0;i<n;i++)a[i]/=n;
        }
    }
    int n,N,m,kind,f[maxn][30],h[30],F[maxn],g[2][maxn];
    complex<double>a[maxn],b[maxn],c[maxn],d[maxn],v[maxn];
    ll ans1[maxn],ans2[maxn],ans3[maxn],ans4[maxn];//
    void multply(complex<double>*x,complex<double>*y,ll *z){
        for(int i=0;i<N;i++)v[i]=x[i]*y[i];
        fft::idft(v,N);
        for(int i=0;i<n;i++)z[i]=(ll)(v[i].real()+0.5);
    }
    void MTT(int *x,int *y,int *z){
        for(int i=0;i<N;i++)a[i]=b[i]=c[i]=d[i]=complex<double>(0,0);//
        for(int i=0;i<n;i++)a[i].real(x[i]>>15);
        for(int i=0;i<n;i++)b[i].real(x[i]&32767);
        for(int i=0;i<n;i++)c[i].real(y[i]>>15);
        for(int i=0;i<n;i++)d[i].real(y[i]&32767);
        fft::dft(a,N);fft::dft(b,N);fft::dft(c,N);fft::dft(d,N);
        multply(a,c,ans1);multply(a,d,ans2);
        multply(b,c,ans3);
        multply(b,d,ans4);
        for(int i=0;i<n;i++)z[i]=(ans1[i]*32768%MOD*32768%MOD+(ans2[i]+ans3[i])*32768%MOD+ans4[i])%MOD;
    }
    int main(){
        n=read();m=read();kind=read();
        f[0][0]=1;h[0]=1;
        int mx=0;
        for(mx=1;h[mx-1]<n;mx++)h[mx]=h[mx-1]*m;mx--;
        for(int i=1;i<=n;i++){
            for(int j=0;j<=mx;j++)if(i-h[j]>=0){
                for(int k=0;k<=j;k++)f[i][j]=(f[i][j]+f[i-h[j]][k])%MOD;
            }
        }
        for(int i=0;i<=n;i++)for(int j=0;j<=mx;j++)F[i]=(F[i]+f[i][j])%MOD;
        int x=0;
        for(int i=0;i<=n;i++)g[x][i]=F[i];
        n++;
        N=1;// 
        while(N<n+n)N*=2;fft::init(N);
        for(int k=1;k<kind;k++)MTT(g[x],F,g[1-x]),x=1-x;
        int sum=0;
        for(int i=0;i<=n;i++)sum=(sum+g[x][i])%MOD;
        printf("%d",sum);
        return 0;
    }
    View Code

    注意清空complex的时候实部和虚部(imag)要一起清空。

    好看的模板:L_0_Forever_LF

    【快速数论变换】NTT

    只记录重要概念,证明略去。

    (一)原根

    当(a,m)=1时,对于满足a^x=1(%m)的最小正整数x,称x为a模m的阶。

    根据欧拉定理a^φ(m)=1(%m),当x=φ(m)时,称a为m的原根。

    以下只讨论m为素数的情况,则当a为m的原根时,a^0~a^(p-2)取遍1~p-1所有值。

    模m有原根的充要条件:m=2,4,p^e,2*p^3,p是奇素数。(也就是说,m为素数时一定有原根)

    求m的原根:p-1= p1^a1 * p2^a2 * pk^ak,g是p的原根当且仅当对于所有的pi满足g^[ (p-1)/pi ] ≠ 1 (%p)

    例题:【51NOD】1135 原根

    (二)快速数论变换

    当模数为形如p=r*2^k+1的素数(费马素数)时,则有n|p-1,可以进行NTT。

    先找到p的原根g(p=998244353 || 1004535809,g=3)

    在原来FFT的基础上,omega[i]=g^[ (p-1)/n*i ] % p,倒数为逆元。

    IDFT时,除以n改为乘n的逆元。

    (三)模数任意的NTT

    找到三个费马素数满足相乘结果>n*(m-1)^2,分别进行NTT后用CRT合并。

    p=998244353,1004535809,469762049,g=3。

    我写的是FFT合并。

    (四)离散对数

    当(a,p)=1时,若满足a^x=b (%p),则称在模p意义下,x是b以a为底的离散对数,即logab=x(单个快速求解可用BSGS算法)。

    1.对于x*y=z(%p),有log x+log y=log z(%p-1),因此离散对数常用于乘法转加法(生成函数)。

    2.对于x^y=z(%p),有y*log x=log z(%p-1)。

    其中log以p的原根g为底。

    例题:【BZOJ】3992: [SDOI2015]序列统计 NTT+生成函数

    【生成函数】母函数

    生成函数的三大要素:①选择项,②大小,③元素个数。一般最终要求某个“大小”的元素个数。

    对于一类组合对象构成的集合A:

    1.每个元素a∈A都定义了一个非负整数的”大小“,记为|a|。

    2.大小为n的元素个数记为$A_n$。

    那么A的一般生成函数是

    $$A(x)=sum_{i=0}^{n}A_ix^i$$

    在这里每个元素都抽象成”大小“,元素a可以理解为有|a|的单位元素的元素。

    组合对象集合D为A和B的笛卡尔积,即D中的每个元素都是A中某个元素a和B中某个元素b组成的有序二元组(a,b),那么显然有D(x)=A(x)B(x)。

    若干一般生成函数的乘积中,第n项的含义是:每个选择项取一个元素,大小相加为n的元素个数。

    每个生成函数本质上是一个集合,那么若干生成函数的乘积就是★每个集合取一个元素的组合。例如生成函数A,B,C,A*B*C的每个元素就是有序三元组(a,b,c)。

    指数型生成函数是

    $$A(x)=sum_{i=0}^{n}A_ifrac{x^i}{i!}$$

    这样就会有:

    $$D_n=sum_{i+j=n}A_iB_jfrac{(i+j)!}{i!j!}=sum_{i+j=n}A_iB_jinom{i+j}{i}$$

    这里乘(i+j)!是因为这只是系数,后面要除以(i+j)!。

    理解为每个元素内部有序就可以了,这样元素内部是排列。

    生成函数都是处理对于n个选择项各选一个组成对应”大小“的元素个数,而一般生成函数元素内部是组合,指数型生成函数元素内部是排列

    一般生成函数还有个化简公式,令x∈[-1,1]时套等比数列公式即可收敛:

    $$sum_{i=0}^{n}x^i=frac{1}{1-x}$$

    指数型生成函数也有个化简公式——泰勒展开:

    $$sum_{i=0}^{n}frac{x^i}{i!}=e^x$$

    例题:

    1.热身:苹果只能取偶数个,橘子只能取1~4个,求拿n个水果的方案数

    题解:定义每种水果为一个集合(每个集合选一个),“大小”为水果个数,最后求总集合”大小“为n的元素个数。

    f(x)=(1+x^2+x^4+...)*(1+x+x^2+x^3+x^4)。

    如果觉得集合很难理解,不妨用”选择项“这个词。

    2.【BZOJ】3771: Triple FTT+生成函数

    题意:给定n个物品,价值为ai,物品价格互不相同,求选一个或两个或三个的价值为x的方案数,输出所有存在的x和对应方案数。ai<=40000。

    题解:要求什么就定义什么为”大小“,所以定义”大小“为价值,[第一个物品][第二个物品][第三个物品]为三个选择项。

    那么每个选择项的每个系数记录对应价值的物品数量(1个)。

    这样拼起来就好了吗?不是,物品不能重复取,所以拼起来之后再容斥掉选相同的。

    我们可以直接写出选一个物品的集合的生成函数f,两个相同物品的g和三个相同物品的h。

    考虑有AAB,ABA,BAA,AAA四种不合法情况,答案就是f^3-3f*g+2h。最后这个求得排列数,需要/3!。选1个或2个的随便推推也一样。

    3.【BZOJ】3992: [SDOI2015]序列统计 NTT+生成函数

    题意:给定一个[0,m-1]范围内的数字集合S,从中选择n个数字(可重复)构成序列。给定x,求序列所有数字乘积%m后为x的序列方案数%1004535809。1<=n<=10^9,3<=m<=8000,m为素数,1<=x<=m-1

    题解:要求乘积,定义”大小“为数字的乘积。但是我们不能加减”大小“啊?

    换成离散对数就可以了,然后定义每个数字为选择项,答案就是f^n。

    4. [母函数]HDU 1521——排列组合 

    题意:有n种物品,并且知道每种物品的数量。要求从中选出m件物品的排列数。例如有两种物品A,B,并且数量都是1,从中选2件物品,则排列有”AB”,”BA”两种。

    题解:定义”大小“为物品数量,选择项为每种物品,那么组合数就是一般生成函数(元素内部有序,严格按物品编号排),排列数就是指数型生成函数(元素内部带标号,可以打乱)。

    这里指数型生成函数,指的是方案就是元素内部带标号的方案数,这样再计算过程中那个公式自动/i!后进行计算再乘回(i+j)!的。

    暴力枚举求解这个生成函数。

    另一部分知识:生成函数公式的化简

    参考:什么是生成函数? by M67

    首先根据二项式定理:

    $$(1+x)^n=sum_{i=0}^{n}inom{n}{i}x^i$$

    扩展到负数和实数,即广义二项式定理,只要将组合数表示成下降幂即可:

    $$(1+x)^{n}=sum_{i=0}^{infty}frac{n^{underline{i}}}{i!}x^i$$

    这里有一个特殊的变换,当n>0时:

    $$(1-x)^{-n}=sum_{i=0}^{infty}frac{(-n)^{underline{i}}}{i!}(-x)^i$$

    $$(1-x)^{-n}=sum_{i=0}^{infty}frac{(-1)^i*n^{overline{i}}}{i!}*(-1)^i*x^i$$

    $$(1-x)^{-n}=sum_{i=0}^{infty}inom{n+i-1}{n-1}*x^i$$

    所以这是将i个相同的数分割n个非空部分的方案数的生成函数。

    另外还常用等比数列递推公式来收敛等比数列:

    $$sum_{i=0}^{infty}(x^q)^i=frac{1}{1-x^q}$$

    另外对于有限的等比数列还常用类似错位相减的方法,左边乘上(1-x)就会变成右边分子,即:

    $$sum_{i=0}^{n}x^i=frac{1-x^{n+1}}{1-x}$$

    例题:

    1.求:

    $$g(x)=(1+x^2+x^4+...)(1+x^5+x^{10}+...)(1+x+x^2+x^3+x^4)(1+x)$$

    用上面提到的技巧即可得到:

    $$g(x)=frac{1}{1-x^2}*frac{1}{1-x^5}*frac{1-x^5}{1-x}*(1+x)$$

    不断地约分,最后用平方差公式化简可得:

    $$g(x)=frac{1}{(1-x)^2}=(1-x)^{-2}=sum_{i=0}^{infty}(i+1)x^i$$ 

    再代入上面那个结论就可以得到多项式了。

    2.食物:用上面的技巧化简约分最后剩个小组合数,n=10^500用读入取模。

    【多项式求逆】

    核心原理是$\%x^{frac{n}{2}}$的多项式平方后可以转化为$\%x^n$。

    已知多项式f(x),求多项式g(x),满足f(x)g(x)=1(%x^n)。

    $$f(x)g(x)=1(\%x^{frac{n}{2}})$$

    $$f(x)h(x)=1(\%x^n)$$

    现在已知g(x),求h(x)。

    $$f(x)g(x)-1=0(\%x^{frac{n}{2}})$$

    将1移到左边后就可以平方了。

    $$f(x)^2g(x)^2-2f(x)g(x)+1=0(\%x^n)$$

    将1换成f(x)*g(x),从而将h(x)代入,提出f(x)消去。

    $$f(x)g(x)^2-2g(x)+h(x)=0(\%x^n)$$

    最终得到:

    $$h(x)=g(x)(2-f(x)g(x))(\%x^n)$$

    然后就可以递归求解,边界条件:当n=1时,f(x)g(x)=1(%x),g(0)是f(0)在%MOD意义下的数论逆元。

    注意:每次n不同都要重新预处理Omega[]。                                                                                                                              

    例题:【BZOJ】4555: [Tjoi2016&Heoi2016]求和 排列组合+多项式求逆 或 斯特林数+NTT

    【拉格朗日插值法】

    参考:拉格朗日插值法(图文详解) by Angel_Kitty

    对于一个n次多项式,如果已知n+1个点,可以构造拉格朗日多项式L(x):

    $$L(x)=sum_{i=0}^{n}y_il_i(x)$$

    其中$l_i(x)$为插值基函数:

    $$l_i(x)=prod_{j=0,j eq i}^{n}frac{x-x_j}{x_i-x_j}$$

    通过代入需要的x即可得到答案。

    每次插值的复杂度为O(n^2)。

    例题:【BZOJ】4559: [JLoi2016]成绩比较 计数DP+排列组合+拉格朗日插值

    当横坐标连续时,上式可以表示为:

    $$f(x)=sum_{i=0}^{n}f(i)*prod_{j=0,j eq i}^{n}frac{x-j}{i-j}$$

    预处理阶乘的逆元$frac{1}{i!}$,预处理$v=prod x-j$,每次乘上v/(i-j),分母是两段阶乘,再根据负数的个数判断正负性。

    复杂度O(n)。

    例题:【BZOJ】3453: tyvj 1858 XLkxc 拉格朗日插值(自然数幂和)

  • 相关阅读:
    UITableView学习笔记
    IOS基础之设置APP的名字、设置图标、添加等待加载时的图片
    UIScrollView,UIPageControl
    UIPickerView基本用法
    最大公约数和最小公倍数
    快速幂、快速乘
    素数筛
    最小生成树
    BZOJ1070 [SCOI2007]修车
    BZOJ1109 [POI2007]堆积木Klo
  • 原文地址:https://www.cnblogs.com/onioncyc/p/8410969.html
Copyright © 2011-2022 走看看