代码改变世界

P3803 【模板】多项式乘法(FFT)

2019-04-22 17:52  一只弱鸡丶  阅读(285)  评论(0编辑  收藏  举报

传送门:

参考博客 1:大佬  attack

参考博客 2:大佬  胡小兔

在这里再膜拜一下这两位大佬 Orz%%%

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
#define re register
#define pb push_back
#define mp make_pair
#define P pair<int,int>
const int N=1e7+10;
const double Pi=acos(-1.0);
void read(int &a)
{
    a=0;
    int d=1;
    char ch;
    while(ch=getchar(),ch>'9'||ch<'0')
        if(ch=='-')
            d=-1;
    a=ch-'0';
    while(ch=getchar(),ch>='0'&&ch<='9')
        a=a*10+ch-'0';
    a*=d;
}
void write(int x)
{
    if(x<0)
        putchar(45),x=-x;
    if(x>9)
        write(x/10);
    putchar(x%10+'0');
}
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
}a[N],b[N];
int r[N],limit=1,l=0;
complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
inline void FFT(complex *A,int type)
{
    for(re int i=0;i<limit;i++)
        if(i<r[i])
            swap(A[i],A[r[i]]);
    for(re int mid=1;mid<limit;mid<<=1)
    {
        complex Wn(cos(Pi/mid),type*sin(Pi/mid));///mid只是当前要做的区间的一半
        for(int R=mid<<1,j=0;j<limit;j+=R)///以间隔为R 为分块
        {
            complex w(1,0);
            for(re int k=0;k<mid;k++,w=Wn*w)///枚举左区间,同时搞出右区间
            {
                complex x=A[j+k],y=A[j+k+mid];
                A[j+k]=x+w*y;
                A[j+k+mid]=x-w*y;
            }
        }
    }
}
int main()
{
    int n,m;
    read(n);
    read(m);
    for(re int i=0;i<=n;i++)
    {
        int x;
        read(x);
        a[i].x=x;
    }
    for(re int i=0;i<=m;i++)
    {
        int x;
        read(x);
        b[i].x=x;
    }
    while(limit<=n+m)
        limit<<=1,l++;
    for(re int i=0;i<limit;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    FFT(a,1);
    FFT(b,1);
    for(re int i=0;i<limit;i++)
        a[i]=a[i]*b[i];
    FFT(a,-1);
    for(re int i=0;i<=n+m;i++)
        printf("%d ",(int)(a[i].x/limit+0.5));
    return 0;
}