Link
设(n=65536,x=omega_n^k)。
下面是利用倒推法得到的大致过程:
( ext{FWT})得到(a_i=1),( ext{CU})得到(a_i=x^i),( ext{IDFT})得到(a_i=[i=k]),( ext{QR})得到答案。
( ext{FWT})
( ext{FWT})的过程是从小到大对(iin [0,16))做( ext{CR(}i,i,M ext{)}),其中(M=egin{pmatrix}1&1\1&-1end{pmatrix})。
为了保证(M)是酉矩阵,令(M=egin{pmatrix}frac1{sqrt 2}&frac1{sqrt 2}\frac1{sqrt 2}&-frac1{sqrt 2}end{pmatrix})即可。
( ext{IDFT})
首先是( ext{reverse}),或者说用(overline{omega_n^i})替代(omega_n^i)。如果我们忽视该操作,那么最后调用( ext{QR})返回的就是(n-k),并不会影响我们求出答案。当然直接用(overline{omega_n^i})做也不会有问题。
然后是( ext{bitreverse}),只需将( ext{CU(}i,2^i ext{)})改成( ext{CU(}i,2^{15-i} ext{)})即可。
接着是( ext{DFT})过程,从小到大枚举(iin [0,16)),然后先对所有满足(2^isubseteq S)进行(a_S=a_Somega_{2^{i+1}}^S)操作,再调用( ext{CR(}i,i,M ext{)})。
第一个操作可以对(jin[0,i))做( ext{CR(}i,j,N ext{)}),其中(N=egin{pmatrix}1&0\0&omega_{2^{i+1}}^{2^j}end{pmatrix})。
最后是每个位置除以(n),忽略该操作不会造成影响。
Tips
在( ext{FWT})和( ext{IDFT})中,我们使用的(M)矩阵相对于正确的(M)矩阵是乘上了一个(frac1{sqrt2})的系数。
这样正好可以达成IDFT最后一步的效果。
#include<cmath>
#include<numeric>
#include"unnamable.h"
using cp=std::complex<double>;
const int N=65536;const double pi=acos(-1),s2=sqrt(0.5);
cp omega(int k){return cp(cos(2*pi*k/N),sin(2*pi*k/N));}
cp FWT[2][2]={{s2,s2},{s2,-s2}},FFT[2][2]={{1,0},{0,0}};
cp SOL(int t)
{
INI(16);
for(int i=0;i<16;++i) CR(i,i,FWT);
for(int i=0;i<16;++i) CU(i,1<<(15-i));
for(int i=0;i<16;CR(i,i,FWT),++i) for(int j=0;j<i;++j) FFT[1][1]=omega(1<<(15-i+j)),CR(i,j,FFT);
return omega(-QR());
}