[FFT]快速傅里叶变换（分治版）

最新推荐文章于 2024-07-31 00:30:00 发布

Agakiss

最新推荐文章于 2024-07-31 00:30:00 发布

阅读量1.1k

点赞数 2

分类专栏：学习笔记文章标签： FFT 快速傅里叶变换分治FFT

本文链接：https://blog.csdn.net/qq_15647239/article/details/90677933

版权

学习笔记专栏收录该内容

15 篇文章 1 订阅

订阅专栏

FFT真好恰！

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<complex>
using namespace std;
#define cp complex<double>
#define ll long long
#define PI acos(-1.0)
#define MAXN 4000010

cp a[MAXN], b[MAXN], c[MAXN];
int n, m, lim;

inline ll read() {
    ll s = 0, w = 1;
    char c = getchar();
    for (; !isdigit(c); c = getchar()) if (c == '-') w = -1;
    for (; isdigit(c); c = getchar()) s = (s << 1) + (s << 3) + (c ^ 48);
    return s * w;
}

cp omega(int n, int k) {
    return cp{cos(2 * PI * k / n), sin(2 * PI * k / n)};
}
void fft(cp *a, int n, bool inv) {
    if (n == 1) return;
    static cp buf[MAXN];
    int m = n / 2;
    for (register int i = 0; i < m; i++) {
        buf[i] = a[2 * i];
        buf[i + m] = a[2 * i + 1];
    }
    for (register int i = 0; i < n; i++)
        a[i] = buf[i];
    fft(a, m, inv);
    fft(a + m, m, inv);
    for (register int i = 0; i < m; i++) {
        cp x = omega(n, i);
        if (inv) x = conj(x);
        buf[i] = a[i] + x * a[i + m];
        buf[i + m] = a[i] - x * a[i + m];
    }
    for (register int i = 0; i < n; i++)
        a[i] = buf[i];
}
int main() {
    n = read(), m = read();
    for (register int i = 0; i <= n; i++)
        a[i] = {read(), 0};
    for (register int i = 0; i <= m; i++)
        b[i] = {read(), 0};
    int lim = 1;
    while (lim <= n + m) lim *= 2;
    for (int i = n + 1; i <= lim; i++) a[i] = {0, 0};
    for (int i = m + 1; i <= lim; i++) b[i] = {0, 0};
    fft(a, lim, true), fft(b, lim, true);
    for (register int i = 0; i <= lim; i++)
        c[i] = a[i] * b[i];
    fft(c, lim, false);
    for (register int i = 0; i <= n + m; i++)
        printf("%d ", (int)((c[i].real() / lim) + 0.5));
    return 0;
}