搬运一下之前在luogu blog写的。
参考http://picks.logdown.com/posts/177631-fast-fourier-transform
废话不写了。首先补零将多项式项数补为2的整数次幂。下标均从0开始。
DFT（递归版）
将系数表示转为点值表示。设项数为n，即求多项式将所有1的n次根代入所得的值。
设多项式为A。ωn为n次根。
A(ωn^m)=A0(ω(n/2)^m)+ωn^mA1(ω(n/2)^m)
A(ωn^(m+n/2))=A0(ω(n/2)^m)-ωn^mA1(ω(n/2)^m)
A0为将A的偶数项取出组成一个新的n/2-1次多项式，A1为奇数项同理。
到只有一项的时候显然将1代入所得值就是那个系数。
一个分治与合并的过程。这样可以递归实现。
DFT（迭代版）
递归常数太大空间爆炸。考虑递归的过程。
0 1 2 3 4 5 6 7
0 2 4 6 | 1 3 5 7
0 4 | 2 6 | 1 5 | 3 7
0 | 4 | 2 | 6 | 1 | 5 | 3 | 7
000 100 010 110 001 101 011 111
可以发现二进制下反过来就是从小到大排序的了。
于是开始按照这种顺序排好序就可以愉快的合并了。
每次迭代出来的东西就是其代表的那个多项式的点值表示。
细节：如何求二进制下翻转并排序
r[i]表示i的翻转。有r[i]=(r[i>>1]>>1)|(i&1)*(n>>1)
即考虑它的翻转与右移一位时有何不同。r[i>>1]>>1即其去掉最后一位的翻转。然后给它添上第一位，i&1即这位是1还是0，n>>1即其最高位代表的值。
排序时直接扫一遍，若i<r[i]交换两个位置的数。正确性因为显然有r[r[i]]=i。
模板：
struct complex{ double x,y; complex operator +(const complex&a) const { return (complex){x+a.x,y+a.y}; } complex operator -(const complex&a) const { return (complex){x-a.x,y-a.y}; } complex operator *(const complex&a) const { return (complex){x*a.x-y*a.y,x*a.y+y*a.x}; } };//复数类 void DFT(int n,complex *a,int p) { for (int i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]); for (int i=2;i<=n;i<<=1) { complex wn=(complex){cos(2*PI/i),p*sin(2*PI/i)};//i项式用i次根 p为1 DFT 为-1 IDFT for (int j=0;j<n;j+=i) { complex w=(complex){1,0}; for (int k=j;k<j+(i>>1);k++,w=w*wn) { complex x=a[k],y=w*a[k+(i>>1)]; a[k]=x+y,a[k+(i>>1)]=x-y; } } } } for (int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1)*(n>>1);
IDFT
将点值表示化为系数表示。其实并没有太搞懂不过毕竟写的时候只要传个参数那就先放放吧。注意最后值除以项数。
NTT
没啥区别就是模意义下的。那用原根代替复数就好了。
原根随便都能求出来。
模数需要为2^n*k+1的形式且为质数。因为需要求得2^i次根，也即原根的(p-1)/2^i次，这个东西显然需要是整数。
常用模数有
998244353=2^23*119+1
1004535809=2^21*479+1
469762049=2^26*7+1
都能跑几百万项。
IDFT时用原根的逆元。最后乘项数的逆元。
模板
for (int i=2;i<=n;i<<=1) { int wn=pow(o,(p-1)/i); for (int j=0;j<n;j+=i) { int w=1; for (int k=j;k<j+(i>>1);k++,w=1ll*w*wn%p) { int x=a[k],y=1ll*w*a[k+(i>>1)]%p; a[k]=(x+y)%p,a[k+(i>>1)]=(x-y+p)%p; } } }