FFT

Posted 2021-03-17 ryedii-blog

多项式

[ egin{aligned} A(x)& =sum_{i=0}^na_ix_i &=(x_0,y_0),(x_1,y_1),cdots,(x_n,y_n) end{aligned} ]

卷积

[ A(x)B(x)=sum_{i=0}^nsum_{j=0}^na_ib_jx^{i+j} ]

若 (Aast B=C)，可知求 (C) 的系数数列时间复杂度 (O(n^2))。

考虑点值表达式的卷积

[ A(x)=(x_0,y_0),(x_1,y_1),cdots,(x_n,y_n) ]

[ B(x)=(x_0,y_0'),(x_1,y_1'),cdots,(x_n,y_n')) ]

可得出 (C) 的一个不完整表达

[C(x)=(x_0,y_0y_0'),(x_1,y_1y_1'),cdots,(x_n,y_ny_n')]

若将点对个数增至 (2n+1) 个，可知求出 (C) 的点值数列时间复杂度 (O(n))。

DFT & IDFT

DFT(Discrete Fourier Transform，离散傅里叶变换) 是系数数列 (a_0,a_1,cdots,a_n) 转为点值数列的过程，IDFT(Inverse Discrete Fourier Transform，逆离散傅里叶变换) 是 DFT 的逆运算。

FFT

考虑 DFT 和 IDFT 时间复杂度难以接受，可以利用 FFT(Fast Fourier Transform，快速傅里叶变换) 将其加速。

单位复根

(omega^n=1) 在复数域下有 (n) 个不同的根。

[ e^{ix}=cos x+isin x e^{2pi i}=1=omega^n \omega=e^{frac{2pi i}n} ]

定义主次单位根

[ omega_n=e^{frac{2pi i}n} ]

对于其它单位根，记作

[ omega_n^k=e^{frac{2pi ik}n} ]

存在以下性质

[ omega_{2n}^{k}+omega_{2n}^{n+k}=0 ]

[ (omega_{2n}^k)^2=omega_{n}^k ]

考虑用单位复根的性质优化 DFT。

对于多项式 (A(x)=sum_{i=0}^{2n-1}a_ix^i)，其中 (2n=2^k)，我们将 (omega_{2n}^0,omega_{2n}^1,cdots,omega_{2n}^{{2n}-1}) 代入公式计算点值。

现在重定义两个多项式

[ A_0(x)=a_0+a_2x+a_4x^2+cdots+a_{2n}x^n ]

[ A_1(x)=a_1+a_3x+a_5x^2+cdots+a_{2n-1}x^n ]

显然

[ A(x)=A_0(x^2)+xA_1(x^2) ]

将单位复根代入上式

[ egin{aligned} A(omega_{2n}^k)&=A_0((omega_{2n}^k)^2)+omega_{2n}^kA_1((omega_{2n}^k)^2)&=A_0(omega_n^k)+omega_{2n}^kA_1(omega_n^k)A(omega_{2n}^{n+k})&=A_0((omega_{2n}^{n+k})^2)+omega_{2n}^{n+k}A_1((omega_{2n}^{n+k})^2)&=A_0(omega_n^k)-omega_{2n}^kA_1(omega_n^k) end{aligned} ]

发现对于 (kin[0,1,cdots,n-1]) (A(omega_{2n}^k)) 和 (A(omega_{2n}^{n+k})) 是可以在一起计算的。

于是有以下伪代码

function FFT(complex A[], int Length)
    if Length == 1 return
    complex A0[Length / 2], A1[Length / 2]
    for int i = 0 to Length / 2 - 1 with i += 1
        A0[i] = A[i * 2]
        A1[i] = A[i * 2 + 1]
    FFT(A0, Length / 2)
    FFT(A1, Length / 2)
    complex wn = (cos(2 * Pi / Length), sin(2 * Pi / Length))
    complex w = (1, 0)
    for int i = 0 to Length / 2 - 1 with i += 1, w *= wn
        A[i] = A0[i] + A1[i] * w
        A[i + Length / 2] = A0[i] - A1[i] * w

考虑 IDFT，IDFT 是 DFT 的逆运算，令 (DFT({a_i})={b_i})，已知

[b_k=sum_{i=0}^{n-1}a_iomega_n^{ki}]

存在结论

[a_k=frac 1 nsum_{i=0}^{n-1}b_iomega_n^{-ki}]

证明：将前式代入后式

[ egin{aligned} a_k&=frac 1 nsum_{i=0}^{n-1}b_iomega_n^{-ki}&=frac 1 nsum_{i=0}^{n-1}omega_n^{-ki}sum_{j=0}^{n-1}a_jomega_n^{ij}&=frac 1 nsum_{j=0}^{n-1}a_jsum_{i=0}^{n-1}omega_n^{-ki}omega_{n}^{ij}&=frac 1 nsum_{j=0}^{n-1}a_jsum_{i=0}^{n-1}omega_{n}^{i(j-k)} end{aligned} ]

考虑 (sum_{i=0}^{n-1}omega_n^{i(j-k)})

当 (j=k)，
[sum_{i=0}^{n-1}omega_n^{i(j-k)}=sum1=n]

当 (j eq k)，

[ egin{aligned} sum_{i=0}^{n-1}omega_n^{i(j-k)}&=sum_{i=0}^{n-1}(omega_n^{j-k})^i&=frac{1-(omega_n^{j-k})^n}{1-omega_n^{j-k}}&=frac{1-(omega_n^n)^{j-k}}{1-omega_n^{j-k}}=frac{1-1}{1-omega_n^{j-k}}=0 end{aligned} ]

所以，

[ egin{aligned} a_k&=frac 1 nsum_{j=0}^{n-1}a_jsum_{i=0}^{n-1}omega_{n}^{i(j-k)}&=frac 1 ncdot na_k=a_k end{aligned} ]

发现 DFT 和 IDFT 的公式形式一样，对于指数上的 (-1) 只需要在代码中加入一个开关即可。

蝴蝶定理

考虑 FFT 的分治过程，以 (n=16) 为例

[{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}]

[{0,2,4,6,8,10,12,14},{1,3,5,7,9,11,13,15}]

[{0,4,8,12},{2,6,10,14},{1,5,9,13},{3,7,11,15}]

[{0,8},{4,12},{2,10},{6,14},{1,9},{5,13},{3,11},{7,15}]

其二进制下表示，

[{0000,1000},{0100,1100},{0010,1010},{0110,1110},]

[{0001,1001},{0101,1101},{0011,1011},{0111,1111}]

观察发现，若干次蝴蝶操作（奇偶分治）后，其数位二进制翻转后是连续的。

对于二进制翻转，可用递推计算，rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0)。

考虑用蝴蝶定理使 FFT 的过程避免递归，可以先将 ({a_i}) 按 (rev) 序重新排列，在分治树上从下往上，

function FFT(complex A[], int Length, int on)
    for int i = 0 to Length - 1 with i += 1
        if i < Rev[i]
            swap(A[i], A[Rev[i]])
    for int k = 1 to n - 1 with k *= 2
        complex wn = (cos(Pi / k), sin(on * Pi / k))
        for int i = 0 to n with i += k * 2
            complex w = (1, 0)
            for int j = i to i + k - 1 with j += 1
                complex x = a[j]
                complex y = a[j + k]
                a[j] = x + y
                a[j + k] = x - y
    if on == -1
        for int i = 0 to n - 1 with i += 1
            a[i] /= n

然后枚举当前合并的长度 (2k)，枚举合并区间开始位置 (i)，枚举区间中的元素 (a_j)。

---

Luogu P3803 多项式乘法(FFT)

#include <iostream>
#include <cmath>

using namespace pl;

const int N = 4200000 + 7;
const double PI = acos(-1);

int n, m;
complex a[N], b[N];
int rev[N];

void mul(complex a[], int n, complex b[], int m);
void fft(complex a[], int n, int on);

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; ++i) scanf("%lf", &a[i].r);
    for (int i = 0; i <= m; ++i) scanf("%lf", &b[i].r);
    mul(a, n, b, m);
    for (int i = 0; i <= n + m; ++i)
        printf("%d", (int)(a[i].r + 0.5));

    return 0;
}

void mul(complex a[], int n, complex b[], int m) {
    for (n = m = n + m - 1; n ^ (n & -n); n += n & -n);
    fft(a, n, 1), fft(b, n, 1);
    for (int i = 0; i < n; ++i) a[i] = a[i] * b[i];
    fft(a, n, -1);
}

void fft(complex a[], int n, int on) {
    static int rev[N], lstn;
    if (n != lstn) {
        lstn = n;
        for (int i = 0; i < n; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1 ? n >> 1 : 0);
    }
    complex t, w, x;
    for (int k = 1; k < n; k <<= 1) {
        t = (complex){ cos(PI / k), on * sin(PI / k) };
        for (int i = 0; i < n; i += k << 1) {
            w = (complex){ 1, 0 }；
            for (int j = i; j < i + k; ++j, w = w * t) 
                x = w * a[j + k], a[j + k] = a[j] - x, a[j] = a[j] + x;
        }
    }
    if (on == -1)
        for (int i = 0; i < n; ++i) a[i].r /= n;
}

---

[mathcal{by Ryedii}]