拆系数FFT学习笔记
拆系数FFT学习笔记
拆系数FFT
当题目中取模的数不是NTT模数的时候，我们无法利用原根来进行快速数论变换，这个时候就要用到毛啸论文里提到的拆系数FFT。
大致思路
拆系数FFT实际上是将多项式卷积之后的值具体算出来，普通的FFT由于精度误差较大，当然无法胜任。
于是可以考虑将一个较大的数拆成两个较小的数，即将\(x\)拆成\(a \times t + b\)，这样相当于将\(x\)给\(t\)进制划分，在常数\(t\)取根号级别的时候，\(a,b\)的大小也是根号级别的，这样我们再用精度较高的浮点数就可以直接FFT来运算了。
例如：
将\(a_1b_1,(a_1b_2+a_2b_1),a_2b_2\)对应算出来再各自乘上相应的系数即可。
这样需要做8次dft，发现中间有两个部分系数相同，对于这两个系数相同的部分我们可以直接点值操作之后再做一次dft即可，这样dft的次数就降到了7次了。
优化
DFT
上面的方法好像常数很大，很多时候需要用到两次DFT合并为一次的方法。
考虑一个多项式A的DFT本质上求的是\(A(\omega_n^k)=\sum_{j=0}^{n-1}A_{j}\times \omega_n^{kj}\)，现在我们需要同时求出A和B，也就是算出\(A(\omega_n^k),B(\omega_n^k)\)。
假设我们知道了\(F_p[k]=A(\omega_n^k)+iB(\omega_n^k),F_q[k]=A(\omega_n^k)-iB(\omega_n^k)\)，我们就可以通过\(F_p[k]\)和\(F_q[k]\)来直接求出\(A(\omega_n^k),B(\omega_n^k)\)了。
于是为了合并两次DFT，我们希望\(F_p\)和\(F_q\)之间能够存在某一种关系，使得求出了某一个之后可以快速求出另外一个。
于是我们来推导\(F_p\)和\(F_q\)的关系，为了方便下面的公式的书写，我们约定\(X=\frac{2\pi kj}{n}\)。
我们发现\(F_p[n-k]\)和\(F_q[k]\)在本质上就是共轭的，于是通过\(F_p\)来计算\(F_q\)就非常方便了。
简单分析一下原因，由于\(\omega_{n}^{n-k}=\omega_n^{-k}\)，而\(\omega_{n}^{k}\)和\(\omega_n^{-k}\)实部和虚部的元素都是相同的，只是符号相反而已，\(F_p[k]\)和\(F_q[k]\)的系数刚好也是符号相反，也就导致\(F_p[k],F_p[n-k],F_q[k],F_q[n-k]\)四个项的它们本质的结果只是一些元素的线性组合，这样从\(F_p\)推出\(F_q\)也就不是那么难了。
IDFT
求出了点值之后我们可以快速将点值乘起来，但是要怎么IDFT回去？
对于两个点值序列\(A,B\)，我们仍然可以用刚才提到的方法，将每一项合并成\(A+iB\)，将新的多项式IDFT之后得到的系数的实部和虚部就对应了真正\(A\)和\(B\)了。
/*=======================================
* Author : ylsoi
* Time : 2019.3.14
* Problem : luogu4245
* E-mail : ylsoi@foxmail.com
* ====================================*/
#include<bits/stdc++.h>
#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define debug(x) cout<<#x<<"="<<x<<" "
#define fi first
#define se second
#define mk make_pair
#define pb push_back
#define double long double
typedef long long ll;
using namespace std;
void File(){
freopen("luogu4245.in","r",stdin);
freopen("luogu4245.out","w",stdout);
}
template<typename T>void read(T &_){
_=0; T fl=1; char _c=getchar();
for(;!isdigit(_c);_c=getchar())if(_c=='-')fl=-1;
for(;isdigit(_c);_c=getchar())_=(_<<1)+(_<<3)+(_c^'0');
_*=fl;
}
const int maxn=1e5+10;
const double pi=acos(-1);
int n,m,mod,a[maxn<<2],b[maxn<<2],c[maxn<<2];
namespace Poly{
struct cp{
double x,y;
cp(double xx=0,double yy=0){x=xx,y=yy;}
cp conj(){return cp(this->x,-this->y);}
friend cp operator + (cp a,cp b){return (cp){a.x+b.x,a.y+b.y};}
friend cp operator - (cp a,cp b){return (cp){a.x-b.x,a.y-b.y};}
friend cp operator * (cp a,cp b){return (cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
friend cp operator * (cp a,double b){return cp(a.x*b,a.y*b);}
friend cp operator / (cp a,double b){return cp(a.x/b,a.y/b);}
}om[maxn<<2],aa[maxn<<2],bb[maxn<<2],cc[maxn<<2],dd[maxn<<2];
int lim,cnt,dn[maxn<<2];
inline void dft(cp *A,int ty){
if(ty==-1)reverse(A+1,A+lim);
REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
for(int len=1;len<lim;len<<=1){
cp w=om[len<<1];
for(int L=0;L<lim;L+=len<<1){
cp wk=cp(1,0);
REP(i,L,L+len-1){
cp u=A[i],v=A[i+len]*wk;
A[i]=u+v;
A[i+len]=u-v;
wk=wk*w;
}
}
}
if(ty==-1)REP(i,0,lim-1)A[i]=A[i]/lim;
}
void mtt(int *A,int *B,int *C,int la,int lb){
REP(i,0,lim-1){
dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
aa[i]= i<=la ? cp(A[i]>>15,A[i]&32767) : cp(0,0);
bb[i]= i<=lb ? cp(B[i]>>15,B[i]&32767) : cp(0,0);
}
dft(aa,1),dft(bb,1);
REP(i,0,lim-1){
int j=(lim-i)&(lim-1);
cp a1=(aa[i]+aa[j].conj())*cp(0.5,0);
cp a0=(aa[i]-aa[j].conj())*cp(0,-0.5);
cp b1=(bb[i]+bb[j].conj())*cp(0.5,0);
cp b0=(bb[i]-bb[j].conj())*cp(0,-0.5);
cc[i]=a1*b1+a1*b0*cp(0,1);
dd[i]=a0*b1+a0*b0*cp(0,1);
}
dft(cc,-1),dft(dd,-1);
REP(i,0,lim-1){
int a1=(ll)(cc[i].x+0.5)%mod;
int b1=(ll)(cc[i].y+0.5)%mod;
int c1=(ll)(dd[i].x+0.5)%mod;
int d1=(ll)(dd[i].y+0.5)%mod;
C[i]=(((ll)a1<<30)+((ll)(b1+c1)<<15)+(ll)d1)%mod;
}
}
}
int main(){
File();
read(n),read(m),read(mod);
REP(i,0,n)read(a[i]);
REP(i,0,m)read(b[i]);
using namespace Poly;
lim=1,cnt=0;
while(lim<=n+m)lim<<=1,++cnt;
if(!cnt)cnt=1;
REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
for(int i=lim;i;i>>=1)
om[i]=cp(cos(2.0*pi/i),sin(2.0*pi/i));
Poly::mtt(a,b,c,n,m);
REP(i,0,n+m)printf("%d ",(c[i]+mod)%mod);
return 0;
}
