【FFT优化】[UOJ#34]多项式乘法

题目描述

这是一道模板题。

给你两个多项式，请输出乘起来后的多项式。

输入格式

第一行两个整数 n 和 m，分别表示两个多项式的次数。

第二行 n+1 个整数，分别表示第一个多项式的 0 到 n 次项前的系数。

第三行 m+1 个整数，分别表示第一个多项式的 0 到 m 次项前的系数。

输出格式

一行 n+m+1 个整数，分别表示乘起来后的多项式的 0 到 n+m 次项前的系数。

样例一

input

1 21 21 2 1

output

1 4 5 2

explanation

(1+2x)⋅(1+2x+x2)=1+4x+5x2+2x3。

限制与约定

0≤n,m≤105，保证输入中的系数大于等于 0 且小于等于 9。

时间限制：1s

空间限制：256MB

分析

本博客主要用于存FFT优化模板.
对于两个实系数的多项式A(x),B(x)，长度均为N(N为2的整次数幂)，令

P(x)=A(x)+iB(x)Q(x)=A(x)−iB(x)

Fp[k],Fq[k]为P(x),Q(x)进行DFT后的值.
令X=2πjkN

Fp[k]=A(ωkN)+iB(ωkN)=∑j=0N−1AjωjkN+i∑j=0N−1BjωjkN=∑j=0N−1(Aj+iBj)∗(cos(X)+isin(X))

Fq[k]=A(ωkN)−iB(ωkN)=∑j=0N−1AjωjkN−i∑j=0N−1BjωjkN=∑j=0N−1(Aj−iBj)∗(cos(X)+isin(X))=∑j=0N−1Aj∗cos(X)+iAjsin(X)−iBjcos(X)+Bjsin(X)=∑j=0N−1(Ajcos(X)+Bjsin(X))+i(Ajsin(X)−Bjcos(X))=∑j=0N−1conj((Ajcos(X)+Bjsin(X))−i(Ajsin(X)−Bjcos(X)))=∑j=0N−1conj((Ajcos(−X)−Bjsin(−X))+i(Ajsin(−X)+Bjcos(−X)))=∑j=0N−1conj(Aj(cos(−X)+isin(−X))+iBj(cos(−X)+isin(−X)))=∑j=0N−1conj((Aj+iBj)∗(cos(−X)+isin(−X))=∑j=0N−1conj(Ajω−jkN+iBj(ω−jkN))=conj(A(ω−kN)+iB(ω−kN))=conj(A(ωN−kN)+iB((ωN−kN))=conj(Fp[N−k])

这样我们就可以通过一次DFT求出Fp,Fq

DFT(A[k])=Fp[k]+Fq[k]2DFT(B[k])=Fp[k]−Fq[k]2i=−iFp[k]−Fq[k]2

如果只是做多项式乘法C(x)=A(x)∗B(x)，直接使N变为A,B长度的和，然后

DFT(C[k])=(Fp[k]2−Fq[k]2)∗−i4=(Fp[k]2−Fp[N−k]2)∗−i4

代码

#include<cstdio>#include<algorithm>#include<cmath>#include<cstring>#define MAXN 300000using namespace std;const double pi=acos(-1);int n,m,N,ans[MAXN+10];char s[MAXN+10];struct cpx{    double r,i;    inline cpx(){    }    inline cpx(double r,double i):r(r),i(i){    }    inline cpx operator+(const cpx &b)const{        return cpx(r+b.r,i+b.i);    }    inline cpx operator-(const cpx &b)const{        return cpx(r-b.r,i-b.i);    }    inline cpx operator*(const cpx &b)const{        return cpx(r*b.r-i*b.i,r*b.i+i*b.r);    }    inline cpx &operator*=(const cpx &b){        return *this=*this*b;    }    inline cpx conj()const{        return cpx(r,-i);    }}a[MAXN+10],b[MAXN+10],t[MAXN+10];void fft(cpx *a,int N,int f){    int i,j=0,k,t;    for(i=1;i<N-1;i++){        for(t=N;j^=t>>=1,~j&t;);        if(i<j)            swap(a[i],a[j]);    }    for(i=1;i<N;i<<=1){        cpx wn(cos(pi/i),f*sin(pi/i));        t=i<<1;        for(j=0;j<N;j+=t){            cpx w(1,0);            for(k=0;k<i;k++,w*=wn){                cpx x(a[j+k]),y(w*a[j+k+i]);                a[j+k]=x+y,a[j+k+i]=x-y;            }        }    }    if(f==-1)        for(i=0;i<N;i++)            a[i].r/=N;}void mul(cpx *a,cpx *b,cpx *c,int N){    int i,j;    for(i=0;i<N;i++)        t[i]=cpx(a[i].r,b[i].r);    fft(t,N,1);    for(i=0;i<N;i++){        j=(N-i)&(N-1);        c[i]=(t[i]*t[i]-(t[j]*t[j]).conj())*cpx(0,-0.25);    }    fft(c,N,-1);    // c refer t    //a[i]=(c[i]+c[k-i])/2    //b[i]=-i(c[i]-c[k-i])/2    //a[i]*b[i]=(c[i]*c[i]-c[k-i]*c[k-i])*((-i)/4）}void read(){    int lenb,lena,i;    scanf("%s",s);    lena=strlen(s);    for(i=0;i<lena;i++)        a[i].r=s[lena-i-1]-'0';    scanf("%s",s);    lenb=strlen(s);    for(i=0;i<lenb;i++)        b[i].r=s[lenb-i-1]-'0';    for(N=1;N<=lena+lenb-2;N<<=1);}void print(){    int i;    for(i=0;i<N;i++){        ans[i]+=a[i].r+0.2;        ans[i+1]+=ans[i]/10;        ans[i]%=10;    }    for(i=N-1;i;i--)        if(ans[i])            break;    for(;i>=0;i--)        printf("%d",ans[i]);    puts("");}int main(){    read();    mul(a,b,a,N);    print();}