搬运一下之前在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
     
      >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
      
       >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
     
      >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;            }        }    }