算法学习笔记：FFT方法解析

算法学习笔记

FFT(FastFourierTransformation)" role="presentation">FFT(FastFourierTransformation)$FFT(FastFourierTransformation)$，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
~~我默认~~大家都知道多项式乘法O(n2)" role="presentation">O(n2)$O(n^2)$的方法，FFT" role="presentation">FFT$FFT$就是把它变成O(nlogn)" role="presentation">O(nlogn)$O(nlogn)$。

首先，FFT" role="presentation">FFT$FFT$要先将系数表示法转化为点值表示法。

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

接着，FFT" role="presentation">FFT$FFT$要用数学方法转化两种表示法。

然后，FFT" role="presentation">FFT$FFT$要先DFT" role="presentation">DFT$DFT$然后再逆DFT" role="presentation">DFT$DFT$。

离散傅里叶变换（DiscreteFourierTransform" role="presentation">DiscreteFourierTransform$DiscreteFourierTransform$，缩写为DFT" role="presentation">DFT$DFT$），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT" role="presentation">DTFT$DTFT$的频域采样。在实际应用中通常采用快速傅里叶变换计算DFT" role="presentation">DFT$DFT$。
假设f(x)=a0+a1x+a2x2+…+an−1xn−1" role="presentation">f(x)=a0+a1x+a2x2+…+an?1xn?1$f(x)=a_0+a_{1}x+a_2{x}^2+\dots +a_{n?1}{x}^{n?1}$
构造f0(x)=a0+a2x+a4x2+…+an−2xn2−1" role="presentation">f0(x)=a0+a2x+a4x2+…+an?2xn2?1$f_0(x)=a_0+a_{2}x+a_4{x}^2+\dots +a_{n?2}{x}^{\frac{n}{2}?1}$
构造f1(x)=a1+a3x+a5x2+…+an−1xn2−1" role="presentation">f1(x)=a1+a3x+a5x2+…+an?1xn2?1$f_1(x)=a_1+a_{3}x+a_5{x}^2+\dots +a_{n?1}{x}^{\frac{n}{2}?1}$
那么f(x)=f0(x2)+xf1(x2)" role="presentation">f(x)=f0(x2)+xf1(x2)$f(x)=f_0(x^2)+x{f}_1(x^2)$
假设xk=ωnk" role="presentation">xk=ωkn$x_k={\omega }_n^k$，那么xk+n2=−xk" role="presentation">xk+n2=?xk$x_{k+\frac{n}{2}}=?x_k$
这样，我们就可以用f0(x)" role="presentation">f0(x)$f_0(x)$和f1(x)" role="presentation">f1(x)$f_1(x)$的点值表示法求出f(x)" role="presentation">f(x)$f(x)$的点值表示法
时间O(nlogn)" role="presentation">O(nlogn)$O(nlogn)$。
至于逆DFT" role="presentation">DFT$DFT$，水证一下。
对于∑ωipn" role="presentation">∑ωipn$\frac{\sum {\omega }^{ip}}{n}$，
当p=0" role="presentation">p=0$p=0$时，原式=1" role="presentation">=1$=1$。
当p≠0" role="presentation">p≠0$p\ne 0$时，原式=1n1−qn1−q=0" role="presentation">=1n1?qn1?q=0$=\frac{1}{n}\frac{1?q^n}{1?q}=0$。
所以，原式=[p=0]" role="presentation">=[p=0]$=[p=0]$。
用最原始的式子：
(17)(18)ci=∑j∑kajbk[−i+j+k=0](19)=1n∑j∑kajbk∑lω−ilωjlωkl(20)=1n∑lω−il(∑jajωjl)(∑kbkωkl)" role="presentation">ci=∑j∑kajbk[?i+j+k=0]=1n∑j∑kajbk∑lω?ilωjlωkl=1n∑lω?il(∑jajωjl)(∑kbkωkl)(17)(18)(19)(20)$$\begin{array}{}\text{(17)}& \text{(18)}& c_i& =\sum _j\sum _k{a}_j{b}_k[?i+j+k=0]\text{(19)}& & =\frac{1}{n}\sum _j\sum _k{a}_j{b}_k\sum _l{\omega }^{?il}{\omega }^{jl}{\omega }^{kl}\text{(20)}& & =\frac{1}{n}\sum _l{\omega }^{?il}(\sum _j{a}_j{\omega }^{jl})(\sum _k{b}_k{\omega }^{kl})\end{array}$$
也就是说先将a" role="presentation">a$a$，b" role="presentation">b$b$各做一次DFT" role="presentation">DFT$DFT$，再将a" role="presentation">a$a$、b" role="presentation">b$b$相乘得到d" role="presentation">d$d$，最后将d" role="presentation">d$d$做一次逆DFT" role="presentation">DFT$DFT$。

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

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

算法学习结束。

题目

2179:FFT快速傅立叶" role="presentation">2179:FFT快速傅立叶$2179:FFT快速傅立叶$

TimeLimit:10Sec" role="presentation">TimeLimit:10Sec$TimeLimit:10Sec$
MemoryLimit:259MB" role="presentation">MemoryLimit:259MB$MemoryLimit:259MB$

Description" role="presentation">Description$Description$

给出两个n" role="presentation">n$n$位10" role="presentation">10$10$进制整数x" role="presentation">x$x$和y" role="presentation">y$y$，你需要计算x⋅y" role="presentation">x?y$x?y$。

Input" role="presentation">Input$Input$

第一行一个正整数n" role="presentation">n$n$。 第二行描述一个位数为n" role="presentation">n$n$的正整数x" role="presentation">x$x$。 第三行描述一个位数为n" role="presentation">n$n$的正整数y" role="presentation">y$y$。

Output" role="presentation">Output$Output$

输出一行，即x⋅y" role="presentation">x?y$x?y$的结果。

SampleInput" role="presentation">SampleInput$SampleInput$

1
3
4

SampleOutput" role="presentation">SampleOutput$SampleOutput$

12

数据范围：

n≤60000" role="presentation">n≤60000$n\le 60000$

题解

这是FFT" role="presentation">FFT$FFT$的入门题目，只要将大数an−1an−2…a1a0¯" role="presentation">an?1an?2…a1a0ˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉˉ$\overline{a_{n?1}{a}_{n?2}\dots a_1{a}_0}$转化为f(x)=an−1xn−1+an−2xn−2+…+a1x+a0" role="presentation">f(x)=an?1xn?1+an?2xn?2+…+a1x+a0$f(x)=a_{n?1}{x}^{n?1}+a_{n?2}{x}^{n?2}+\dots +a_{1}x+a_0$，然后将两个多项式相乘，并且将x=10" role="presentation">x=10$x=10$代入即可。
时间O(nlogn)" role="presentation">O(nlogn)$O(nlogn)$，空间O(n)" role="presentation">O(n)$O(n)$。

总结

当我们需要求解ci=∑ajbi−j" role="presentation">ci=∑ajbi?j$c_i=\sum a_j{b}_{i?j}$这样的式子的时候，可以利用FFT" role="presentation">FFT$FFT$小优化。当然，如果不用能过，那在好不过了。

标程

#include 
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;
}