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

题目描述

这是一道模板题。

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

输入格式

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

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

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

输出格式

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

样例一

input

1 2
1 2
1 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 300000
using 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();
}