前言
随便乱写,意识流警告,反正只是当个板子
因为要写long double所以可能比普通7次的FFT还要慢?!
站在常数的底端.jpg
MTT
假设我们要求
A
(
x
)
∗
B
(
x
)
A(x)*B(x)
A(x)∗B(x)
先把P拆成
P
\sqrt P
P,每个数拆成
a
∗
p
+
b
a*\sqrt p+b
a∗p+b,然后分别求FFT,最后组合起来
这样要7次
我们可构造
P
(
x
)
=
A
(
x
)
+
i
B
(
x
)
,
Q
(
x
)
=
A
(
x
)
−
i
B
(
x
)
P(x)=A(x)+iB(x),Q(x)=A(x)-iB(x)
P(x)=A(x)+iB(x),Q(x)=A(x)−iB(x)
可以发现
Q
(
ω
k
)
=
c
o
n
j
(
P
(
ω
n
−
k
)
)
Q(\omega^k)=conj(P(\omega^{n-k}))
Q(ωk)=conj(P(ωn−k))
于是我们只需要对P求一遍DFT就知道了对Q求DFT的结果
然后知道了P和Q的点值我们就能还原A和B的点值
然后构造
M
=
A
(
x
)
+
i
B
(
x
)
M=A(x)+iB(x)
M=A(x)+iB(x),IDFT回去可以发现M的实部是A(x)的插值,虚部是B(x)的插值
然后就只用4次了
Code
const db pi=acos(-1);
struct Com{
db r,i;
Com (db _r=0,db _i=0) {r=_r;i=_i;}
friend Com operator + (Com x,Com y) {return Com(x.r+y.r,x.i+y.i);}
friend Com operator - (Com x,Com y) {return Com(x.r-y.r,x.i-y.i);}
friend Com operator * (Com x,Com y) {return Com(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
}t1[N],t2[N],t3[N],t4[N];
Com conj(Com a) {return Com(a.r,-a.i);}
void DFT(Com *a,int len,int flag) {
fo(i,0,len-1) if (i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<len;i<<=1) {
Com wn=Com(cos(pi/i),flag*sin(pi/i));
for(int j=0;j<len;j+=i<<1) {
Com w=Com(1,0);
for(int k=0;k<i;++k,w=w*wn) {
Com u=a[j+k],v=w*a[j+k+i];
a[j+k]=u+v;a[j+k+i]=u-v;
}
}
}
if (flag==-1) fo(i,0,len-1) a[i].r/=len,a[i].i/=len;
}
void mult(int *a,int *b,int *c,int n) {
int lg,len;
for(lg=0,len=1;len<=n;len<<=1) lg++;
fo(i,0,len-1) rev[i]=rev[i>>1]>>1|((i&1)<<(lg-1));
fo(i,0,len-1) t1[i]=t2[i]=t3[i]=t4[i]=Com(0,0);
fo(i,0,n) t1[i]=Com(a[i]&M,a[i]>>15),t2[i]=Com(b[i]&M,b[i]>>15);
DFT(t1,len,1);DFT(t2,len,1);
fo(i,0,len-1) {
int j=(len-i)&(len-1);
Com da=(t1[i]+conj(t1[j]))*Com(0.5,0);// a%r
Com db=(t1[i]-conj(t1[j]))*Com(0,-0.5);// a/r
Com dc=(t2[i]+conj(t2[j]))*Com(0.5,0);// b%r
Com dd=(t2[i]-conj(t2[j]))*Com(0,-0.5);// b/r
t3[i]=da*dc+da*dd*Com(0,1);
t4[i]=db*dc+db*dd*Com(0,1);
}
DFT(t3,len,-1);DFT(t4,len,-1);
fo(i,0,n) {
int ka=(ll)(t3[i].r+0.5)%Mo,kb=(ll)(t3[i].i+0.5)%Mo;
int kc=(ll)(t4[i].r+0.5)%Mo,kd=(ll)(t4[i].i+0.5)%Mo;
c[i]=((ll)ka+((ll)(kb+kc)<<15)+((ll)kd<<30))%Mo;
}
}