快速傅里叶变换

2023/10/12 多项式卷积傅里叶变换

学了数字信号处理里的FFT于是上洛谷看了许久的cpp

思路：

对于朴素的 $O(n^2)$ 的乘法求多项式显然过不了 $1e6$ 的高时间复杂度， $此时引入离散傅里叶变换(DFT)$ 和 $逆离散傅里叶变换$。
我们可以通过简单模拟知道多项式的乘法就是其系数的卷积运算，然后时间域的卷积等于其频域的乘法，即 $f[t] * g[t] = F(jw) .* G(jw) (这里区分 .* 和 * 分别为普通乘法和卷积运算)$。这里的等式是连续函数的，离散条件下该等式也成立。
这里我们过程就是
$A和B的多项式系数→(DFT)→A和B的频域形式相乘→(IDFT)→ANSWER的多项式系数$

这个过程中，我们 $DFT$ 采用 $FFT$ 进行正转换和逆转换，其中的复杂度为 $O(nlogn)$, 频域相乘的复杂度为 $O(n)$。

实现

$FFT$ 实际上就是将一整个序列平均分成两个序列进行 $DFT$，发现这俩个部分存在重复的值，便可以不断实现 $O(logn)$ 的划分，最终实现递归。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2e6 + 10;
const double pi = acos(-1.0);
int n, m, r[N], lim = 1, l;
struct com{  //结构体实现虚数坐标轴
    double x, y;
    com(double xx = 0, double yy = 0){
        x = xx; y = yy;
    }
}a[N], b[N];
// 虚数加减乘法的实现
com operator + (com a, com b) {return com(a.x + b.x, a.y + b.y);}
com operator - (com a, com b) {return com(a.x - b.x, a.y - b.y);}
com operator * (com a, com b) {return com(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void fft(com *x, int type){

    // 由蝶形运算的下标位置对应关系
    for(int i = 0; i < lim; ++ i){
        if(r[i] > i) swap(x[i], x[r[i]]);
    }

    // 枚举每次折中下标
    for(int mid = 1; mid < lim; mid <<= 1){
        com k(cos(pi / mid), type * sin(pi  / mid)) ;
        // 幂次方
        for(int R = mid << 1, j = 0; j < lim; j += R){
            com w(1, 0); 
            // 幂次方的基底
            for(int i = 0; i < mid; ++ i, w = w * k){
                com aa = x[j + i], bb = w * x[j + i + mid];

                // 反向递推
                x[i + j] = aa + bb;
                x[i + j + mid] = aa - bb;
            }
        }
    }
}
int main(){
    cin >> n >> m;
    for(int i = 0; i <= n; ++ i) cin >> a[i].x;
    for(int i = 0; i <= m; ++ i) cin >> b[i].x;
    
    //由基2的FFT原理，我们构造一个最小2的幂次
    while(lim <= n + m) lim <<= 1, l ++;
    

    // 由蝶形运算法，我们求出每个下标序列二进制的倒序
    // 该序列位置的时间域可以求出其二进制倒序位置的频域元素值，IDFT同理也可求出
    for(int i = 0; i < lim; ++ i) 
        r[i] = (r[i >> 1] >> 1 | (i & 1) << (l - 1));

    fft(a, 1);  // 求出 A 多项式的DFT
    fft(b, 1);  // 求出 B 多项式的DFT
    
    // 求出 A B 多项式在频域的相乘
    for(int i = 0; i <= lim; ++ i) a[i] = a[i] * b[i];

    // 求出 A 多项式的IDFT
    fft(a, -1);

    for(int i = 0; i <= n + m; ++ i) {
        int xx = a[i].x / lim + 0.5;
        cout << xx << ' ';
    }
}

LOADING

快速傅里叶变换

思路：

实现