FFT学习笔记

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