模拟考某题一开始由于校内OJ太慢直接拆系数FFT跑不过
后来被神仙婊了一顿之后发现复杂度写炸了改了改随便过
模版题：任意模数NTT
三模数NTT
常数巨大，跑的极慢
拆系数FFT
原理是对于两个多项式$ P=\sum\limits_{i=0}^{n-1}P_ix^i \ \ Q=\sum\limits_{i=0}^{m-1}Q_ix^i$
直接$ FFT$计算会发现值域达到$ 10^{23}$会炸精度
设
$ A=\sum\limits_{i=0}^{n-1}(P_i>>15)x^i \ \ B=\sum\limits_{i=0}^{n-1}(P_i\&32767)x^i$
$ C=\sum\limits_{i=0}^{m-1}(Q_i>>15)x^i \ \ D=\sum\limits_{i=0}^{m-1}(Q_i\&32767)x^i$
我们只要求$ (A*C)<<30,(B*C+A*D)<<15,B*D$这三项的和即可
设一次$ DFT/IDFT$为一次操作
暴力实现需要进行$ 8$次操作
精度问题
如果用$ k$次乘法计算$ w_n^k$会损失大量精度
需要利用三角函数预处理单位根，这样可以用$ double$代替$ long \ double$
优化
参考myy的2016年集训队论文
合并$DFT$
设我们要计算$ DFT_A$和$DFT_B$
令$$ P(x)=A(x)+iB(x) \\ Q(x)=A(x)-iB(x)$$
我们只要计算一次$ DFT_P$就可以推出$ DFT_Q$
推导请参考CMXRYNP'S Blog
有
$DFT_A[i]=\frac{DFT_P[i]+DFT_Q[i]}{2}$
$DFT_B[i]=\frac{DFT_P[i]-DFT_Q[i]}{2i}$
合并$IDFT$
同理
但这里甚至不需要求$ IDFT_Q$
事实上$IDFT_P$的实部和虚部分别对应$ IDFT_A$和$IDFT_B$
这样就把$ 8$次操作减少到$4$次了
代码
#include<ctime>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<vector>
#define l putchar('\n')
#define file(x)freopen(x".in","r",stdin);freopen(x".out","w",stdout)
#define block 32768
#define rt register int
#define ll long long
using namespace std;
inline ll read(){
ll x=;char zf=;char ch=getchar();
while(ch!='-'&&!isdigit(ch))ch=getchar();
if(ch=='-')zf=-,ch=getchar();
while(isdigit(ch))x=x*+ch-'',ch=getchar();return x*zf;
}
void write(ll y){if(y<)putchar('-'),y=-y;if(y>)write(y/);putchar(y%+);}
void writeln(const ll y){write(y);putchar('\n');}
int k,m,n,x,y,z,cnt,ans,p;
namespace any_module_NTT{
vector<int>R;
const double PI=acos(-1.0);
struct cp{
double x,y;
cp operator +(const cp s)const{return {x+s.x,y+s.y};}
cp operator -(const cp s)const{return {x-s.x,y-s.y};}
cp operator *(const cp s)const{return {x*s.x-y*s.y,x*s.y+y*s.x};}
}w[][<<];
void FFT(const int n,vector<cp>&A){
A.resize(n);
for(rt i=;i<n;i++)if(i>R[i])swap(A[i],A[R[i]]);
for(rt i=,s=;i<n;i<<=,s++){
for(rt j=;j<n;j+=i<<){
for(rt k=;k<i;k++){
const register cp x=A[j+k],y=w[s][k]*A[i+j+k];
A[j+k]=x+y,A[i+j+k]=x-y;
}
}
}
}
vector<int>Mul(vector<int>&x,vector<int>&y){ int sz=x.size()+y.size()-,lim=;
while(lim<=sz)lim<<=;R.resize(lim);
for(rt i=;(<<i)<lim;i++)
for(rt j=;j<(<<i);j++)w[i][j]={cos(PI*j/(<<i)),sin(PI*j/(<<i))};
vector<cp>AB(lim),CD(lim),AC(lim),BC(lim);
for(rt i=;i<lim;i++)R[i]=(R[i>>]>>)|(i&)*(lim>>);
for(rt i=;i<x.size();i++)AB[i].x=((ll)x[i])&,AB[i].y=((ll)x[i])>>;
for(rt i=;i<y.size();i++)CD[i].x=((ll)y[i])&,CD[i].y=((ll)y[i])>>;
FFT(lim,AB);FFT(lim,CD);
for(rt i=;i<lim;i++){
static cp na,nb,nc,nd;const int pl=(lim-)&(lim-i);
na=AB[i]+(cp){AB[pl].x,-AB[pl].y},nb=AB[i]-(cp){AB[pl].x,-AB[pl].y};
nc=CD[i]+(cp){CD[pl].x,-CD[pl].y},nd=CD[i]-(cp){CD[pl].x,-CD[pl].y};
const cp v1={0.5,},v2={,-0.5};
na=na*v1;nb=nb*v2;nc=nc*v1;nd=nd*v2;
AC[pl]=na*nc+na*nd*(cp){,};
BC[pl]=nb*nc+nb*nd*(cp){,};
}
FFT(lim,AC);FFT(lim,BC);
vector<int>ans(sz);
for(rt i=;i<sz;i++){
ll v1=AC[i].x/lim+0.5,v2=AC[i].y/lim+BC[i].x/lim+0.5,v3=BC[i].y/lim+0.5;
ans[i]=(ll)((v3%p<<)+(v2%p<<)+v1)%p;
}
return ans;
}
}
using namespace any_module_NTT;
vector<int>A,B;
int main(){
n=read();A.resize(n+);m=read();B.resize(m+);p=read();
for(rt i=;i<=n;i++)A[i]=read();for(rt i=;i<=m;i++)B[i]=read();
A=Mul(A,B);
for(rt i=;i<=n+m;i++)write((A[i]+p)%p),putchar(' ');
return ;
}