快速傅里叶变换（FFT）

最新推荐文章于 2025-05-21 23:13:20 发布

✎﹏ℳ๓无情❦

最新推荐文章于 2025-05-21 23:13:20 发布

阅读量194

点赞数

分类专栏： C++ 算法 # 数学文章标签： c++ 算法开发语言

本文链接：https://blog.youkuaiyun.com/qq_61540355/article/details/126810714

版权

算法同时被 3 个专栏收录

17 篇文章

订阅专栏

C++

16 篇文章

订阅专栏

数学

1 篇文章

订阅专栏

文章目录

暂时只能理解这是干嘛的，以及该怎么用。。。

卷积
给定向量： $a=(a_0,a_1,...,a_{n−1})$ ， $b=(b_0,b_1,...,b_{n−1})$ ；
向量和： $a+b=(a_0+b_0,a_1+b_1,...,a_{n−1}+b_{n−1})$
数量积（内积、点积）： $a⋅b=a_0b_0+a_1b_1+...+a_{n−1}b_{n−1}$
卷积： $a⊗b=(c_0,c_1,...,c_{2n−2})$ ，其中 $c_k=∑_{i+j=k}(a_ib_j)$
卷积的最典型的应用就是多项式乘法（多项式乘法就是求卷积）。

$a⊗b=IDFT_{2n}(DFT_{2n}(a)⋅DFT_{2n}(b))$ ，即： $a⊗b=DFT^{−1}_{2n}(DFT_{2n}(a)⋅DFT_{2n}(b))$

#include <bits/stdc++.h>

#define fp(i, a, b) for (int i = (a), i##_ = (b) + 1; i < i##_; ++i)
using namespace std;
using ll = int64_t;
using db = double;
/*---------------------------------------------------------------------------*/
struct cp
{
    db x, y;
    cp(db real = 0, db imag = 0) : x(real), y(imag){};
    cp operator+(cp b) const
    {
        return {x + b.x, y + b.y};
    }
    cp operator-(cp b) const
    {
        return {x - b.x, y - b.y};
    }
    cp operator*(cp b) const
    {
        return {x * b.x - y * b.y, x * b.y + y * b.x};
    }
};
using vcp = vector<cp>;
using Poly = vector<ll>;
namespace FFT
{
const db pi = acos(-1);
vcp Omega(int L)
{
    vcp w(L);
    w[1] = 1;
    for (int i = 2; i < L; i <<= 1)
    {
        auto w0 = w.begin() + i / 2, w1 = w.begin() + i;
        cp wn(cos(pi / i), sin(pi / i));
        for (int j = 0; j < i; j += 2)
            w1[j] = w0[j >> 1], w1[j + 1] = w1[j] * wn;
    }
    return w;
}
auto W = Omega(1 << 21); // NOLINT
void DIF(cp *a, int n)
{
    cp x, y;
    for (int k = n >> 1; k; k >>= 1)
        for (int i = 0; i < n; i += k << 1)
            for (int j = 0; j < k; ++j)
                x = a[i + j], y = a[i + j + k], a[i + j + k] = (a[i + j] - y) * W[k + j], a[i + j] = x + y;
}
void IDIT(cp *a, int n)
{
    cp x, y;
    for (int k = 1; k < n; k <<= 1)
        for (int i = 0; i < n; i += k << 1)
            for (int j = 0; j < k; ++j)
                x = a[i + j], y = a[i + j + k] * W[k + j], a[i + j + k] = x - y, a[i + j] = x + y;
    const db Inv = 1. / n;
    fp(i, 0, n - 1) a[i].x *= Inv, a[i].y *= Inv;
    reverse(a + 1, a + n);
}
} // namespace FFT

namespace Polynomial
{
// basic operator
void DFT(vcp &a)
{
    FFT::DIF(a.data(), a.size());
}
void IDFT(vcp &a)
{
    FFT::IDIT(a.data(), a.size());
}
int norm(int n)
{
    return 1 << (__lg(n - 1) + 1);
}

// Poly mul
vcp &dot(vcp &a, vcp &b)
{
    fp(i, 0, a.size() - 1) a[i] = a[i] * b[i];
    return a;
}
Poly operator+(Poly a, Poly b)
{
    int maxlen = max(a.size(), b.size());
    Poly ans(maxlen + 1);
    a.resize(maxlen + 1), b.resize(maxlen + 1);
    for (int i = 0; i < maxlen; i++)
        ans[i] = a[i] + b[i];
    return ans;
}
Poly operator*(ll k, Poly a)
{
    Poly ans;
    for (auto i : a)
        ans.push_back(k * i);
    return ans;
}
Poly operator*(Poly a, Poly b)
{
    int n = a.size() + b.size() - 1;
    vcp c(norm(n));
    fp(i, 0, a.size() - 1) c[i].x = a[i];
    fp(i, 0, b.size() - 1) c[i].y = b[i];
    DFT(c), dot(c, c), IDFT(c), a.resize(n);
    fp(i, 0, n - 1) a[i] = int(c[i].y * .5 + .5);
    return a;
}
} // namespace Polynomial
/*---------------------------------------------------------------------------*/
using namespace Polynomial;
const int N = 5e4 + 10;
int main()
{
    ios::sync_with_stdio(0), cin.tie(0);
    int n, m;
    cin >> n >> m;
    Poly a(n + 1), b(m + 1), c(n + m + 1);
    for (int i = 0; i <= n; i++)
        cin >> a[i];
    for (int i = 0; i <= m; i++)
        cin >> b[i];
    c = a * b;
    for (int i = 0; i <= n + m; i++)
        cout << c[i] << " ";
    return 0;
}