算法学习笔记：FFT

最新推荐文章于 2024-10-17 17:45:57 发布

lzx_2014

最新推荐文章于 2024-10-17 17:45:57 发布

阅读量339

点赞数

分类专栏：算法题解 FFT

本文链接：https://blog.csdn.net/lzx_2014/article/details/79690247

版权

算法题解同时被 2 个专栏收录

19 篇文章 0 订阅

订阅专栏

FFT

2 篇文章 0 订阅

订阅专栏

算法学习笔记

$FFT(Fast Fourier Transformation)$ ，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
~~说实话我根本就没看懂，那我就草草地说说。~~
~~我默认~~大家都知道多项式乘法 $O(n^2)$ 的方法， $FFT$ 就是把它变成 $O(nlogn)$ 。

首先， $FFT$ 要先将系数表示法转化为点值表示法。

系数表示法：已知 $n$ 项 $a_i$ ， $f(x)=\sum{a_i x^i}$
点值表示法：已知 $n$ 项 $x_i$ 和 $y_i$ ， $y_i=f(x_i)$
其中， $f(x)$ 为 $n$ 项 $n-1$ 次多项式
可得， $f(x)$ 唯一确定，且两种表示法可相互转化。

接着， $FFT$ 要用数学方法转化两种表示法。

$e^{xi} = cos x + i sin x$
设 $\omega_n = e^{\frac{2\pi i}{n}}$ ，则 $\omega_n^n = (e^{\frac{2\pi i}{n}})^n = e^{2\pi i} = 1$ ， $(\omega_n^k)^n = (\omega_n^n)^k = 1^k = 1$
$(\omega_n^{k+\frac{n}{2}})^2 = \omega_n^{2k+n} = \omega_n^{2k} = (\omega_n^k)^2$

然后， $FFT$ 要先 $DFT$ 然后再逆 $DFT$ 。

离散傅里叶变换（ $Discrete Fourier Transform$ ，缩写为 $DFT$ ），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其 $DTFT$ 的频域采样。在实际应用中通常采用快速傅里叶变换计算 $DFT$ 。
~~说实话我还是没看懂，那我还是草草地说说。~~
假设 $f(x) = a_0 + a_1 x + a_2 x^2 + … + a_{n-1} x^{n-1}$
构造 $f_0(x) = a_0 + a_2 x + a_4 x^2 + … + a_{n-2} x^{\frac{n}{2}-1}$
构造 $f_1(x) = a_1 + a_3 x + a_5 x^2 + … + a_{n-1} x^{\frac{n}{2}-1}$
那么 $f(x) = f_0(x^2) + x f_1(x^2)$
假设 $x_k = \omega_n^k$ ，那么 $x_{k+\frac{n}{2}} = -x_k$
这样，我们就可以用 $f_0(x)$ 和 $f_1(x)$ 的点值表示法求出 $f(x)$ 的点值表示法
时间 $O(nlogn)$ 。
至于逆 $DFT$ ，水证一下。
对于 $\frac{\sum{\omega^{ip}}}{n}$ ，
当 $p=0$ 时，原式 $=1$ 。
当 $p\neq 0$ 时，原式 $=\frac{1}{n}\frac{1-q^n}{1-q}=0$ 。
所以，原式 $=[p=0]$ 。
用最原始的式子：

c i = \sum j \sum k a j b k [- i + j + k = 0] = 1 n \sum j \sum k a j b k \sum l ω - i l ω j l ω k l = 1 n \sum l ω - i l (\sum j a j ω j l) (\sum k b k ω k l) (17) (18) (19) (20)

$\begin{align} \\ c_i &= \sum_j \sum_k {a_j b_k [-i+j+k=0]} \\ &= \frac{1}{n} \sum_j \sum_k {a_j b_k \sum_l {\omega^{-il} \omega^{jl} \omega^{kl}}} \\ &= \frac{1}{n} \sum_l {\omega^{-il} (\sum_j {a_j \omega^{jl}}) (\sum_k {b_k \omega^{kl}})} \\ \end{align}$
也就是说先将

a a $a$ ，

b

$b$ 各做一次

DFT D F T $DFT$ ，再将

a a $a$ 、

b

$b$ 相乘得到

d d $d$ ，最后将

d

$d$ 做一次逆

DFT D F T $DFT$ 。

分治过程中，需要将系数交错排列。
零层： $0, 1, 2, 3, 4, 5, 6, 7$ 。
一层： $0, 2, 4, 6, 1, 3, 5, 7$ 。
两层： $0, 4, 2, 6, 1, 5, 3, 7$ 。
这被叫做蝴蝶变换。

我们为了从下到上计算，必须预处理出这个序列。
将序列转化为二进制数，可以发现神奇的玩意儿。
两层： $000, 100, 010, 110, 001, 101, 011, 111$
零层： $000, 001, 010, 011, 100, 101, 110, 111$
正好倒序，也就是说，可以直接得出。

算法学习结束。
~~写笔记好累的，好奇大佬们都怎么写的。~~

$题目传送门luogu1919$
$题目传送门bzoj2179$

题目

$2179: FFT快速傅立叶$

$Time Limit: 10 Sec$
$Memory Limit: 259 MB$

$Description$

给出两个 $n$ 位 $10$ 进制整数 $x$ 和 $y$ ，你需要计算 $x\cdot y$ 。

$Input$

第一行一个正整数 $n$ 。第二行描述一个位数为 $n$ 的正整数 $x$ 。第三行描述一个位数为 $n$ 的正整数 $y$ 。

$Output$

输出一行，即 $x\cdot y$ 的结果。

$Sample Input$

1
3
4

$Sample Output$

数据范围：

$n\leq 60000$

题解

这是 $FFT$ 的入门题目，只要将大数 $\overline{a_{n-1} a_{n-2} … a_1 a_0}$ 转化为 $f(x) = a_{n-1} x^{n-1} + a_{n-2} x^{n-2} + … + a_1 x + a_0$ ，然后将两个多项式相乘，并且将 $x = 10$ 代入即可。
时间 $O(nlogn)$ ，空间 $O(n)$ 。

总结

当我们需要求解 $c_i = \sum a_j b_{i-j}$ 这样的式子的时候，可以利用 $FFT$ 小优化。~~当然，如果不用能过，那在好不过了。~~

标程

#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
const int N = 1 << 17;
struct Complex{
    double r, i;
    Complex(double _r = 0, double _i = 0) {r = _r; i = _i;}
    Complex operator +(Complex c) {return Complex(r + c.r, i + c.i);}
    Complex operator -(Complex c) {return Complex(r - c.r, i - c.i);}
    Complex operator *(Complex c) {return Complex(r * c.r - i * c.i, r * c.i + i * c.r);}
    Complex operator /(int c) {return Complex(r / c, i / c);}
}a[N], b[N];
int n, m, l, c[N], d[N];
char s[N];
void fft(Complex *a, int p)
{
    for (int i = 0; i < n; i++) if (i < d[i]) swap(a[i], a[d[i]]);
    for (int i = 1; i < n; i <<= 1)
    {
        Complex u = Complex(cos(PI / i), p * sin(PI / i));
        for (int j = 0; j < n; j += (i << 1))
        {
            Complex v = Complex(1, 0);
            for (int k = 0; k < i; k++, v = v * u)
            {
                Complex x = a[j + k], y = v * a[j + k + i];
                a[j + k] = x + y; a[j + k + i] = x - y;
            }
        }
    }
    if (p == -1) for (int i = 0; i < n; i++) a[i] = a[i] / n;
}
int main()
{
    scanf("%d", &n); n--;
    scanf("%s", s); for (int i = 0; i <= n; i++) a[i] = Complex(s[n - i] - '0', 0);
    scanf("%s", s); for (int i = 0; i <= n; i++) b[i] = Complex(s[n - i] - '0', 0);
    m = n << 1; for (n = 1; n <= m; n <<= 1) l++;
    for (int i = 0; i < n; i++) d[i]=((d[i >> 1] >> 1) | (i & 1) << (l - 1));
    fft(a, 1); fft(b, 1);
    for (int i = 0; i < n; i++) a[i] = a[i] * b[i];
    fft(a, -1);
    for (int i = 0; i <= m; i++) c[i] = (int)(a[i].r + 0.5);
    for (int i = 0; i <= m; i++)
        if (c[i] >= 10)
        {
            c[i + 1] += c[i] / 10;
            c[i] %= 10;
            if (i == m) m++;
        }
    for (int p = 0, i = m; i >= 0; i--) if (p || c[i]) {p = 1; printf("%d", c[i]);}
    return 0;
}