FFT详解

最新推荐文章于 2022-03-11 15:53:29 发布

weixin_30492601

最新推荐文章于 2022-03-11 15:53:29 发布

阅读量506

点赞数

原文链接：http://www.cnblogs.com/LiM-817/p/10887264.html

版权

用于计算一类朴素卷积问题。

笔者学习的是这份博客
内容中可能有很多相同之处，敬请谅解。

现在要计算两个一元\(n\)次多项式\(F(x)\)与\(G(x)\)的乘积，如何计算？
前置知识：多项式的表示方法
一. 系数表示法
对于一个\(n\)次多项式\(F(x)\)，它可以被表示成
\[F(x) = a_nx^n+a_{n-1}x^{n-1}+...+a_1x^1+a_0x^0.\]
更加形式化的来说，它可以表示成
\[F(x) = \sum_{i=0}^{n} a_ix^i.\]
举个例子，2次多项式，其中\(a_0=1,a_1=2,a_2=1\)
那么\(F(x)=1x^2+2x+1.\)这样即可通俗的表示出一个\(n\)次多项式。
二. 点值表示法
~~众所周知，两个点确定一个一次函数，三个点确定一个二次函数。~~
所以，\(n+1\)个点确定一个一元\(n\)次多项式。
所以我们可以通过\(n+1\)个点来表示它。
那么相乘之后的点值如何计算？
比如说两个2次多项式\(F(x)=x^2+2x+1\)（红色）与\(G(x)=3x^2-4x-2\)（蓝色），它们的图像如图所示：

那么观察图中\(x=1\)时的情况。此时\(F(1)=4,G(1)=-3.\)
所以，显然，\(F(1)\times G(1)=-12\).也就是说在\(Z=F*G\)这一多项式内，带入\(1\)，得到的结果是\(-12\).
等等，好像有哪里不对。如果说\(Z=F*G\)的话，那么Z的次数应该是\(2n\).
那\(Z\)需要\(2n+1\)个点来确定。但是原来只需要\(n+1\)个点，咋办？
很简单，在原来的多项式里每个都多加\(n\)个点即可。反正多项式已知。
这样就可以用点值来进行操作。也就是说先转成点值，再一乘，再转回来，就是计算流程。
但是好像还是很慢。那么如何优化呢?
复数部分
复数，即形如\(a+bi\)的数，其中\(i^2=-1.\) \(a\)称为实部，\(bi\)称为虚部。
或者说：在一个数轴上（只有x轴），我们可以表示出任何实数。
那么，多加一维（y轴），也就是类似于平面直角坐标系一样，我们就可以表示出任意一个复数。
所以我们把这个坐标系叫做复平面，其中x轴称为实轴，y轴称为虚轴。
复数运算
复数相加：实部相加，虚部相加，例如
\[(a+bi)+(c+di)=(a+c)+(b+d)i.\]
复数相减：同理。
\[(a+bi)-(c+di)=(a-c)+(b-d)i.\]
复数相乘：像一次多项式一样相乘。 注意\(i^2=-1\).
\[(a+bi)(c+di)=ac+(ad+bc)i-bd=(ac-bd)+(ad+bc)i.\]
复数相除：
相信大家都学过共轭根式。同样的，复数也有共轭。
即：\(a+bi\)的共轭为\(a-bi\)。
这两个复数卡乘在一起一定是个实数。即
\[(a+bi)(a-bi)=a^2-(bi)^2=a^2+b^2.\]
所以再除的时候，将分子分母同乘分母的共轭，就可以将分母有理化。
即
\[\frac{a+bi}{c+di}=\frac{(a+bi)(c-di)}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i.\]
复数逆元：
\[\frac{1}{a+bi}=\frac{a}{a^2+b^2}-\frac{b}{a^2+b^2}i.\]
同样的，复数运算在复平面内也有一定规律可循。
考虑在复平面内的两个复数：（借张图）

表示的是复数\((1+4i)\)与复数\((3+2i)\)相乘所得的结果：\((-5+14i)\)。
设其中\((5,0)\)点为位置\(P\)，则\(\angle{POC}=\angle{BOA}.\)
还有：\(\overline{OB}\times\overline{OC}=\overline{OA}.\)
第二个证明：勾股定理。先把\(i\)消掉。
\(\overline{OB}^2=a^2+b^2.\) \(\overline{OC}^2=c^2+d^2.\)
\(\overline{OB}^2\times \overline{OC}^2=(a^2+b^2)(c^2+d^2)=a^2c^2+a^2d^2+b^2c^2+b^2d^2.\)
\(\overline{OA}^2=(ac-bd)^2+(ad+bc)^2=a^2c^2+a^2d^2+b^2c^2+b^2d^2.\)得证。
我们将复数中，复数向量的长度称为模长，向量与x轴正方向的夹角称为幅角。
根据上面的东西：复数相乘时，模长相乘，幅角相加。
单位根
一个n次的单位根即为方程\(x^n=1\)的复数解。
考虑这样一个图：

其中圆上的所有复数模长都是1，这个圆称为单位圆。
考虑\(|x|\)的取值范围：
如果\(|x|<1\)，那么\(|x^n|<1\)（模长*模长，越乘越短）
如果\(|x|>1\)，那么\(|x^n|>1\)（模长*模长，越乘越长）
所以一定有\(|x|=1\)，即点一定在这个单位圆上。
那么还有一条，就是一个n次的单位根共有n个，并且从(1,0)开始，每次逆时针走\(\frac{1}{n}\)个圆周。分别称为\(\omega_{n}^{0},\omega_{n}^{1},...\omega_{n}^{n-1}.\)（逆时针编号！！！）
其实还有\(k\geq n\)和\(k<0\)的单位根，本质上就模个n就可以了，这里不再多加考虑。
有关单位根也有很多性质：

\(\omega_{n}^{0}=1\)
\(\omega_{n}^{k}=(\omega_{n}^{1})^k\)，很好理解，就是每次转\(\frac{1}{n}\)圆周。
\(\omega_{n}^{j} \times \omega_{n}^{k}=\omega_{n}^{j+k}.\)
证明：很简单，代入性质2即可。
根据3，有：
\[\omega_{n}^{j+k}=\omega_{n}^{j} \times \omega_{n}^{k}.\]
\[(\omega_{n}^{k})^{-1}=\frac{1}{\omega_{n}^{k}}=\frac{\omega_{n}^{0}}{\omega_{n}^{k}}=\omega_{n}^{-k}=\omega_{n}^{n-k}.\] （根据上面模个n的说法）
\[\omega_{pn}^{pk}=\omega_{n}^{k}.\]（将蛋糕分成pn份取pk份=将蛋糕分成n份取k份）
\[\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}.\]（相当于绕了半圈\(\rightarrow\)取了个相反数）
\[(\omega_{n}^{k})^j=(\omega_{n}^{j})^k.\]（大小为k，共有j个=大小为j，共有k个）
...差不多就这些了
---
\(FFT\)流程：\(DFT,IDFT\)
\(DFT\)加速版本：分治
首先，你需要~~控制得当~~将整个n-1次多项式（保证n是2的正整数次幂）按系数奇偶拍成两部分：
\[F(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3+...+a_{n-1}x^{n-1}).\]
再设两个n/2项多项式：
\[FL(x)=a_0+a_2x+...+a_{n-2}x^{\frac{n}{2}-1}\]
\[FR(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}\]
那么有\(F(x)=FL(x^2)+x*FR(x^2).\)（不懂可以自己设一下代入进去）
下面，我们来代入单位根，看看会有什么神奇的事情发生。
比如说我们代一个\(\omega_{n}^{k}(k<\frac{n}{2}):\)
\[F(\omega_{n}^{k})=FL((\omega_{n}^{k})^2)+\omega_{n}^{k}\times FR((\omega_{n}^{k})^2)\]
其中，\((\omega_{n}^{k})^2=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}.\)然后代入：
\[F(\omega_{n}^{k})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).\]
然后呢？有啥用？
不难发现，如果我们知道了\(FL(x),FR(x)\)在x取\(\omega_{\frac{n}{2}}^{0},\omega_{\frac{n}{2}}^{1}...\omega_{\frac{n}{2}}^{\frac{n}{2}-1}\)下的所有值，我们就可以知道\(F(x)\)在x取\(\omega_{n}^{0},\omega_{n}^{1}...\omega_{n}^{\frac{n}{2}-1}\)下的所有值（代回原式）。
但是这还不是全部啊，剩下那个部分怎么办？
没关系，考虑\(k<\frac{n}{2}\)，代入\(\omega_{n}^{k+\frac{n}{2}}\)时的情况。

\[F(\omega_{n}^{k+\frac{n}{2}})=FL((\omega_{n}^{k+\frac{n}{2}})^2)+\omega_{n}^{k+\frac{n}{2}}\times FR((\omega_{n}^{k+\frac{n}{2}})^2).\]

其中，\((\omega_{n}^{k+\frac{n}{2}})^2=\omega_{n}^{2k+n}.\)

故
\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k+n}).\]

\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{n}^{2k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{n}^{2k}).\]
（模性质）

\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k+\frac{n}{2}}\times FR(\omega_{\frac{n}{2}}^{k}).\]
（性质3 延展）

\[F(\omega_{n}^{k+\frac{n}{2}})=FL(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}\times FR(\omega_{\frac{n}{2}}^{k}).\]

（取反性质）

然后观察这个式子和前面式子的区别，发现就是把+号改成了-号。也就是说我们可以通过那组数据来算出当\(k>\frac{n}{2}\)时的情况。这就是加速DFT的原理。
IDFT
至此，我们已经可以将其转化为点值，那么如何还原回来呢？使用IDFT插值。
考虑一个向量\(c_0,c_1...c_{n-1}\)，满足\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i.\]
其中\(y_i\)表示乘出来的点的y坐标。
那么，这个式子也就代表一个以\(y\)为系数的方程在取\(x=\omega_{n}^{-k}\)时的复数解。
我们尝试着来化简一波。

\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i\times y_i\]
\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{i})^j)\]
\[c_k=\sum _{i=0}^{n-1}(\omega_{n}^{-k})^i(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i)\]
\[c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j})^i(\omega_{n}^{-k})^i)\]
\[c_k=\sum _{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_{n}^{j-k})^i)\]
\[c_k=\sum _{j=0}^{n-1}a_j\color{red}{(\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i)}\]
发现后面标红的东西就是一个等比数列，那么我们可以直接计算。即

\[\sum_{i=0}^{n-1}X^i=\frac{X^n-1}{X-1}.\]

我们设这个东西为公式\(Z.\)
于是就有
\[c_k=\sum _{j=0}^{n-1}a_jZ(\omega_{n}^{j-k}).\]

观察函数\(Z\)的相关性质。
赫然发现：我们带进去的是单位根！单位根啥性质来着？\(X^n=1!\)
也就是说当\(X≠1\)，也就是\(k≠0\)时，原式\(=0.\) （分子为0）
那么\(X=1\)时，也就是\(k=0\)时呢？这时候这个求和公式就不管用了，带回原式算，发现由于\(1^i=1\)（i是自然数），那么\(Z(x)=n.\)
所以说，什么情况下带进去的值\(k=0\)呢？显然是\(j-k\)的时候！此时单项\(=na_k\).
否则的话会乘一个0，答案固然也是0.
所以有：\[c_k=na_k\]（飞跃性的一步！！！）
回顾整个式子，发现我们算\(c\)的时候连\(a\)的影子都没有，现在我们竟然推出了一个和\(a\)有关的关系式！太神奇了！这就是\(IDFT\)的原理。
至此，\(\text{FFT}\)的重要部分已经全部介绍完毕，下面进入代码实现与细节处理部分。

代码实现
朴素DFT:

complex < double > omega(const int n , const int k) { // 单位根
    complex < double > ret = {cos(2 * Pi / n * k) , sin(2 * Pi / n * k)};
    if(!inv) return ret;
    else return conj(ret);
}

void DFT(complex < double > a , const int n) {
    if(n == 1) return ;
    static complex < double > buf[MAXN];
    const int m = n / 2;
    for(int i = 0; i < m; i++) {
        buf[i] = a[(i << 1)];
        buf[i + m] = a[(i << 1 | 1)];
    }
    copy(buf , buf + n , a);
    complex < double > *a1 = a , *a2 = a + m;
    DFT(a1 , m); DFT(a2 , m);
    for(int i = 0; i < m; i++) {
        complex < double > w = omega(n , i);
        buf[i] = a1[i] + w * a2[i];
        buf[i + m] = a1[i] - w * a2[i];
    }
    copy(buf , buf + n , a);
}

发现这样的效率显然不行，那么考虑迭代形式。
首先观察分治到最后的二进制位表示：

000 001 010 011 100 101 110 111
 0   1   2   3   4   5   6   7
 0   2   4   6 - 1   3   5   7
 0   4 - 2   6 - 1   5 - 3   7
 0 - 4 - 2 - 6 - 1 - 5 - 3 - 7
000 100 010 110 001 101 011 111

然后...我们惊奇的发现最后的二进制位就是之前的二进制位倒过来！TQL!
接着考虑空间优化。
不难发现，由于DFT的两个式子十分相近，所以我们可以只算一次，同时通过技巧性的操作直接覆盖原数组，而不再新开一个。这个操作被称为蝴蝶操作。
蝴蝶操作的代码：

for(int l = 2; l <= n; l <<= 1) {
            int m = l / 2;
            for(complex < double > *p = a; p != a + n; p += l) {
                for(int i = 0; i < m; i++) {
                    complex < double > t = omega[n / l * i] * p[m + i];
                    p[m + i] = p[i] - t;
                    p[i] += t;
                }
            }
        }

然后穿起来，与以前的FFT模板结合即可。

#include <bits/stdc++.h>
#define dbg(x) cerr << #x " = " << x << "\n"
#define INF 0x3f3f3f3f
/* do not give up ! */
/* try your best! */
/* Read the meaning clearly! */
typedef long long LL;
typedef long double ld;
typedef unsigned long long ULL;

using namespace std;

template < typename T > inline void inp(T& t) {
    char c = getchar(); T x = 1; t = 0; while(!isdigit(c)) {if(c == '-') x = -1; c = getchar();}
    while(isdigit(c)) t = t * 10 + c - '0' , c = getchar();t *= x;
}
template < typename T , typename... Args > inline void inp(T& t , Args&... args) {inp(t); inp(args...);}
template < typename T > inline void outp(T t) {
    if(t < 0) putchar('-') , t = -t; T y = 10 , len = 1;
    while(y <= t) y *= 10 , len++; while(len--) y /= 10 , putchar(t / y + 48) , t %= y;
}
template < typename T > inline T mul(T x , T y , T MOD) {x=x%MOD,y=y%MOD;return ((x*y-(T)(((ld)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;}

const int MAXN = 2097152 + 10;
const double Pi = acos(-1.0);
int ans[MAXN];
struct FastFourierTransform {
    complex < double > omega[MAXN] , omegaInverse[MAXN];

    void init(const int n) {
        for(int i = 0; i < n; i++) {
            omega[i] = complex < double > (cos(2 * Pi / n * i) , sin(2 * Pi / n * i));
            omegaInverse[i] = conj(omega[i]);
        }
    }
    void Transform(complex < double > *a , const int n , const complex < double > *omega) {
        int k = 0;
        while((1 << k) < n) ++k;
        for(int i = 0; i < n; i++) {
            int t = 0;
            for(int j = 0; j < k; j++) if(i & (1 << j)) t |= (1 << (k - j - 1));
            if(i < t) swap(a[i] , a[t]);
        }
        for(int l = 2; l <= n; l <<= 1) {
            int m = l / 2;
            for(complex < double > *p = a; p != a + n; p += l) {
                for(int i = 0; i < m; i++) {
                    complex < double > t = omega[n / l * i] * p[m + i];
                    p[m + i] = p[i] - t;
                    p[i] += t;
                }
            }
        }
    }

    void DFT(complex < double > *a , const int n) {
        Transform(a , n , omega);
    }
    void IDFT(complex < double > *a , const int n) {
        Transform(a , n , omegaInverse);
        for(int i = 0; i < n; i++) a[i] /= n;
    }
}FFT;
int main() {
    int nn , mm , x; inp(nn , mm);
    nn++ , mm++;
    int n = 1; while(n < nn + mm) n <<= 1;
    complex < double > a[MAXN] , b[MAXN];
    for(int i = 0; i < nn; i++) {
        inp(x); a[i].real(x);
    }
    for(int i = 0; i < mm; i++) {
        inp(x); b[i].real(x);
    }
    FFT.init(n);
    FFT.DFT(a , n);
    FFT.DFT(b , n);

    for(int i = 0; i < n; i++) a[i] *= b[i];

    FFT.IDFT(a , n);
    for(int i = 0; i < nn + mm - 1; i++) {
        ans[i] = static_cast < int > (floor(a[i].real() + 0.5));
        cout << ans[i] << " ";
    }
    return 0;
}

至此，该代码已经可以效率很高地通过P3803一题！完结撒花！
FFT例题将会另外放在一个帖子~

转载于:https://www.cnblogs.com/LiM-817/p/10887264.html

weixin_30492601

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
FFT详解

用于计算一类朴素卷积问题。笔者学习的是这份博客内容中可能有很多相同之处，敬请谅解。现在要计算两个一元\(n\)次多项式\(F(x)\)与\(G(x)\)的乘积，如何计算？前置知识：多项式的表示方法一. 系数表示法对于一个\(n\)次多项式\(F(x)\)，它可以被表示成\[F(x) = a_nx^n+a_{n-1}x^{n-1}+...+a_1x^1+a_0x^0.\]更加形式化的...
复制链接

扫一扫