最近闲着无聊研究了下(FFT)的常数优化,大概就是各种(3)次变(2or1.5)次之类的,不过没见过啥题卡这个的吧。
关于(FFT)可以看这里:浅谈FFT&NTT。
关于复数
设(x=a+bi),其中(i)是虚数单位,那么我们用(ar x)表示(x)的共轭复数,即(ar x=a-bi)。
共轭复数有一个这样的性质:
证明展开就好了,这个是下面优化的关键。
设(omega_n)为(n)阶单位根,则(overline{omega _n^{x}}=omega_{n}^{-x})。
idft变dft
设(f(x)=sum_{i=0}^{n-1}a_ix^i),注意到:
也就是说我们如果进行一次std::reverse(a+1,a+n)
,然后dft(a)
,在除以(n),我们就完成了一次(idft)。
多项式乘法优化 1
给出多项式(f(x)=sum_{i=0}^{n}a_ix^i,g(x)=sum_{i=0}^{m}b_ix^i),求其卷积。
这里最开始介绍一种非常简洁的优化方法,构造多项式(h(x)):
那么我们只需要取(h^2(x))的虚部除以(2)就是答案,这只需要做两次(FFT)。
多项式乘法优化 2
这个和上面的关联不大,设(X_i)表示多项式(F(x))(dft)之后的系数,(a_i)表示(dft)之前的系数,设(F(x))为(n)项的多项式,且(n=2^k),注意到:
即:(X_i=overline{X_{n-i}})。
这实质上是因为(F)没有虚部的原因,我们换一个有虚部的多项式试试:
等等,我们发现第一个式子和第三个式子很像,两式相加减可以得到:
注意到等式右边就是(a) (dft)完之后的结果,那么对于多项式(F(x),G(x)),我们可以构造一个函数然后(dft)一次,然后(O(n))得到两个多项式(dft)之后的结果,总共只用了一次(FFT)。
当然这个玩意也可以这样用:假设我们现在想求(dft(F(x))),我们把(F(x))奇偶分类,构造多项式:
然后相当于是(0.5)次(FFT)来完成这个事,设(dft(g(x)))每一项为(X_i),(dft(F(x)))每一项为(Y_i),那么推一下可以得到:
注意这里只有(iin [0,n/2))的值,(Y_{n/2})特殊处理一下,后面的可以通过前面得到。
MTT常数优化
( m MTT)就是拆系数( m FFT),设多项式(s(x),t(x)),我们要算(s(x)t(x)),模数任意。
我们拆系数,设拆完了之后是(s(x)=a(x)+b(x)cdot p,t(x)=c(x)+d(x)cdot p)。
构造(F(x)=a(x)+icdot b(x)),(G(x)=c(x)+icdot d(x))。
那么有:
那么相加减可得(a(x),b(x))的(dft)。
令(h(x)={ m dft}(a(x))cdot { m dft}(G(x))={ m dft}(a(x)cdot G(x))={ m dft}(a(x)c(x)+icdot a(x)d(x)))。
那么我们(idft)一次(h(x))就可以得到(a(x)c(x),a(x)d(x))。
同理可以得到(b(x)c(x),b(x)d(x)),一共(4)次(dft)。
代码长这样:
void mul(int *r,int *s,int *t,int len) {
for(N=1,bit=0;N<len;N<<=1,bit++);
for(int i=1;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
for(int i=0;i<N;i++) g[0][i]=cp(r[i]&all,r[i]>>15),g[1][i]=cp(s[i]&all,s[i]>>15);
fft(g[0]),fft(g[1]);
for(int i=0;i<N;i++) {
int j=(N-i)&(N-1);
g[2][j]=(g[0][i]+conj(g[0][j]))*cp(0.5,0)*g[1][i];
g[3][j]=(g[0][i]-conj(g[0][j]))*cp(0,-0.5)*g[1][i];
}fft(g[2]),fft(g[3]);
for(int i=0;i<N;i++) g[2][i]=g[2][i]/N,g[3][i]=g[3][i]/N;
for(int i=0;i<N;i++) {
ll pp=g[2][i].r+0.5,x=g[2][i].i+0.5,y=g[3][i].r+0.5,z=g[3][i].i+0.5;
t[i]=(pp%p+(((x+y)%p)<<15)+((z%p)<<30))%p;
}
}