洛谷P4245：【模板】MTT （拆系数+FFT）

题目传送门：https://www.luogu.org/problemnew/show/P4245

---

题目分析：为什么同一道题的题解我写了两篇blog？当然是因为这是两种不同的方法骗访问量啊。

今天终于写了传闻已久的拆系数+FFT。大概就是令$sp=\lceil \sqrt{p}\rceil$，然后将两个多项式按除以$sp$和模$sp$的结果分成两个多项式。最后将这四个多项式两两相乘，再加权即可。由于FFT的计算结果可能达到$10^{14}$，要预处理单位复数根以保证精度。

其实上面的算法框架不难，如何降低DFT的次数才是关键。如果对四个多项式分别进行DFT，两两相乘后的结果再IDFT，需要8倍DFT的常数。myy的论文中给出了一种用一次DFT同时计算两个实序列在单位复数根处的点值的方法，可以将4次DFT优化到2次。而只要令$F_P=DFT_A+i*DFT_B$，再对$F_P$做IDFT，$A(x)$和$B(x)$就会分别变成$P(x)$的实部和虚部。利用这个方法就可以将IDFT的次数同样减小为2次。

由于三模数NTT需要做9次模意义下的DFT，而拆系数+FFT只需要做4次复数意义下的DFT，后者的速度完爆前者。本题中我写的代码，用前一种方法运行总耗时约12500ms，而后一种方法只要2200ms。

---

插一句题外话，今晚sc给了我一道很妙很妙的题。大概是说如何用1~12辨别两个1~920的数字x,y。我们发现$C_{12}^6=924$，于是可以将每个数映射到一种选取方案，然后记第一个x选了而y没有选的物品。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
#include<vector>
using namespace std;

const int maxn=1000000;
const double pi=acos(-1.0);
typedef long long LL;

struct Complex
{
    double X,Y;
    Complex (double a=0.0,double b=0.0) : X(a),Y(b) {}
} ;

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);}

Complex P[maxn];

vector <Complex> w[maxn];
int Rev[maxn];
int N,Lg;

Complex A[maxn];
Complex B[maxn];
Complex C[maxn];
Complex D[maxn];

int F[maxn];
int G[maxn];

int n,m;
LL p,sp;

void DFT(Complex *a,double f)
{
    for (int i=0; i<N; i++)
        if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);

    for (int len=2; len<=N; len<<=1)
    {
        int mid=(len>>1);
        for (Complex *p=a; p!=a+N; p+=len)
            for (int i=0; i<mid; i++)
            {
                Complex temp=w[mid][i];
                temp.Y*=f;
                temp=temp*p[mid+i];
                p[mid+i]=p[i]-temp;
                p[i]=p[i]+temp;
            }
    }
}

void FFT(int *a,Complex *b,Complex *c)
{
    for (int i=0; i<N; i++) P[i]=Complex((double)(a[i]/sp),(double)(a[i]%sp));
    DFT(P,1.0);
    for (int i=0; i<N; i++)
    {
        c[i]=P[(N-i)%N];
        c[i].Y*=-1.0;
    }
    for (int i=0; i<N; i++)
    {
        b[i]=P[i]+c[i];
        b[i].X/=2.0; b[i].Y/=2.0;
        c[i]=P[i]-c[i];
        c[i].X/=2.0; c[i].Y/=2.0;
        swap(c[i].X,c[i].Y);
        c[i].Y*=-1.0;
    }
}

int main()
{
    freopen("4245.in","r",stdin);
    freopen("4245.out","w",stdout);

    scanf("%d%d%I64d",&n,&m,&p);
    for (int i=0; i<=n; i++) scanf("%d",&F[i]),F[i]%=p;
    for (int i=0; i<=m; i++) scanf("%d",&G[i]),G[i]%=p;

    N=1,Lg=0;
    while (N<n+m+4) N<<=1,Lg++;
    for (int i=0; i<N; i++)
        for (int j=0; j<Lg; j++)
            if (i&(1<<j)) Rev[i]|=(1<<(Lg-j-1));

    int len=1;
    while ((len<<1)<=N)
    {
        for (int i=0; i<len; i++)
        {
            double ang=(double)i*pi/len;
            w[len].push_back( Complex( cos(ang) , sin(ang) ) );
        }
        len<<=1;
    }

    sp=1;
    while (sp*sp<p) sp++;
    FFT(F,A,B);
    FFT(G,C,D);

    for (int i=0; i<N; i++)
    {
        Complex a=A[i]*C[i],b=A[i]*D[i]+B[i]*C[i],c=B[i]*D[i];
        A[i]=a;
        B[i]=b;
        C[i]=c;
    }

    for (int i=0; i<N; i++)
    {
        B[i].Y*=-1.0;
        swap(B[i].X,B[i].Y);
        P[i]=A[i]+B[i];
    }
    DFT(P,-1.0);
    DFT(C,-1.0);

    for (int i=0; i<=n+m; i++)
    {
        P[i].X/=((double)N);
        P[i].Y/=((double)N);
        C[i].X/=((double)N);
        LL ans=(long long)floor( P[i].X+0.5 )%p*(sp*sp%p);
        ans=(ans+ (long long)floor( P[i].Y+0.5 )%p*sp%p )%p;
        ans=(ans+ (long long)floor( C[i].X+0.5 )%p )%p;
        printf("%I64d ",ans);
    }
    printf("\n");

    return 0;
}