fft相关的复习

Posted 2021-02-12 foreverpiano

任意长度卷积 CZT

就是一波推导
[ egin{aligned} b_i &= sum_{j=0}^{n-1} omega^{ij}a_j &= sum_{j=0}^{n-1} omega^{frac{i^2+j^2-(i-j)^2}{2}}a_j &= omega^{frac{i^2}{2}} sum_{j=0}^{n-1}omega^{frac{-(i-j)^2}{2}} a_j omega^{j^2} end {aligned} ]

后面是一个减法卷积，就可以扩展到2的幂次直接fft就好了。

2次dft计算卷积

考虑有两个长度为(n = 2^k)的序列(a(x), b(x))，我们要计算他们的dft。

构造序列(p_k = a_k + ib_k, ; q_k = a_k - ib_k)，

有结论(dft_q(k) = conj(dft_p((n - k) mod n)))。展开，考虑几何意义？？？

我们可以解出(dft_a, dft_b?)。

再做一遍idft就可以了

拆系数fft

记(M = sqrt {mod}?)，把(x?)表示成(x = a imes M + b, b < M?)。

((a imes M + b)(c imes M + d) = ac imes M^2 + (ad + bc) imes M + bd)

对每一项分开算，做7次dft就可以了。

套用上述介绍做法4次dft就够了。

实现上注意在idft的时候，直接把一个序列放在real，另一个放在imag，idft回来直接/N后计算贡献就好了。

以及我们可以直接在一个for里面做解出AB，reverse序列的事情。

下面是关键部分的代码。

poly realmain(poly a, poly b) {
    int n = a.size(), m = b.size();
    prepare(n + m - 1);
    for (int i = 0; i < n; i++) A[i] = cpx(a[i] & 32767, a[i] >> 15);
    for (int i = 0; i < m; i++) B[i] = cpx(b[i] & 32767, b[i] >> 15);
    dft(A, fft_n); dft(B, fft_n);
    for (int i = 0; i < fft_n; i++) {
        int j = (fft_n - i) % fft_n;
        cpx ax, ay, bx, by;
        ax = (A[i] + A[j].conj()) * cpx(0.5, 0);
        ay = (A[i] - A[j].conj()) * cpx(0, -0.5);
        bx = (B[i] + B[j].conj()) * cpx(0.5, 0);
        by = (B[i] - B[j].conj()) * cpx(0, -0.5);
        C[j] = ax * bx + ay * by * cpx(0, 1.0);
        D[j] = ay * bx + ax * by * cpx(0, 1.0);
    }
    dft(C, fft_n); dft(D, fft_n);
    poly ans(n + m - 1, 0);
    for (int i = 0; i < ans.size(); i++) {
        lo ax = lo(C[i].x / fft_n + 0.5) % mod;
        lo ay = lo(C[i].y / fft_n + 0.5) % mod;
        lo bx = lo(D[i].x / fft_n + 0.5) % mod;
        lo by = lo(D[i].y / fft_n + 0.5) % mod;
        ans[i] = ax + ((by + bx) << 15) + (ay << 30);
        ans[i] = (ans[i] % mod + mod) % mod;
    }
    return ans;
}