题目传送门:https://www.luogu.org/problemnew/show/P4245
题目分析:为什么同一道题的题解我写了两篇blog?当然是因为这是两种不同的方法骗访问量啊。
今天终于写了传闻已久的拆系数+FFT。大概就是令 sp=⌈p–√⌉ s p = ⌈ p ⌉ ,然后将两个多项式按除以 sp s p 和模 sp s p 的结果分成两个多项式。最后将这四个多项式两两相乘,再加权即可。由于FFT的计算结果可能达到 1014 10 14 ,要预处理单位复数根以保证精度。
其实上面的算法框架不难,如何降低DFT的次数才是关键。如果对四个多项式分别进行DFT,两两相乘后的结果再IDFT,需要8倍DFT的常数。myy的论文中给出了一种用一次DFT同时计算两个实序列在单位复数根处的点值的方法,可以将4次DFT优化到2次。而只要令 FP=DFTA+i∗DFTB F P = D F T A + i ∗ D F T B ,再对 FP F P 做IDFT, A(x) A ( x ) 和 B(x) B ( x ) 就会分别变成 P(x) P ( x ) 的实部和虚部。利用这个方法就可以将IDFT的次数同样减小为2次。
由于三模数NTT需要做9次模意义下的DFT,而拆系数+FFT只需要做4次复数意义下的DFT,后者的速度完爆前者。本题中我写的代码,用前一种方法运行总耗时约12500ms,而后一种方法只要2200ms。
插一句题外话,今晚sc给了我一道很妙很妙的题。大概是说如何用1~12辨别两个1~920的数字x,y。我们发现 C612=924 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;
}