快速傅里叶变换(FFT)

快速傅里叶变换:在 $O(nlog(n))$ 内求出两个多项式的卷积

###前置知识

####多项式的点表示法

对于任意一个多项式 $A(x)=a_0+a_1x^1+a_2x^2+a_3x^3+...+a_nx^n$,我们都可用 $n+1$ 个点将它表示出来

*** 证明 ***

任取 $n+1$ 个点 $(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2))...(x_n,A(x_n))$,将它代入多项式可得

$$\begin{cases}a_0+a_1x_0+a_2x_0^2+...+a_nx_0^n=A(x_0)\\a_0+a_1x_1+a_2x_1^2+...+a_nx_1^n=A(x_1)\\...\\a_0+a_1x_n+a_2x_n^2+...+a_nx_n^n=A(x_n)\end{cases}$$

通过 (n+1) 个点,我们可以将这个多项式的系数求出来,且这个解与多项式是一一对应的

复数

定义

对于一个复数,我们写成 $a+bi$ 的形式,其中 $a$ 表示实部, $bi$表示虚部,$i=\sqrt{-1}$

基本运算

对于两个复数 $x=a_1+b_1i,y=a_2+b_2i$,有

$x+y=(a_1+a_2)+(b_1+b_2)i$

$x-y=(a_1+a_2)+(b_1+b_2)i$

$x\times y=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i$

几何意义

对于任意一个复数，我们都可以将它在坐标系中用向量表示

对于相加与相减，他的几何意义与向量相同

对于相乘,得到的复数的模长等于两个复数的模长和,与 $x$ 轴的夹角等于两个复数的夹角和

复数的单位根

把一个单位圆分成 $n$ 份

对于单位根上的每一份,我们记作 $w_n^k$,它的角度就为 $\frac{2k\pi}{n}$,模长为1

性质

$\forall i\ne j, w_n^i\ne w_n^j$

$w_n^0=w_n^n$

$w_n^k=cos(\frac{2k\pi}{n})+sin(\frac{2k\pi}{n})i$

$w_{2n}^{2k}=w_n^k$

$w_n^{k+\frac{n}{2}}=-w_n^k$

证明过程

对于两个多项式 $A(x)={a_0+a_1x+...+a_nx^n},B(x)={b_0+b_1x+...+b_mx^m}$

令 $C(x)=A(x)B(x)$

很明显, $C(x)$ 一共有 $n+m+1$ 项,所以我们只需要找到 $n+m+1$ 个点,就可以用点表示法表示出 $C(x)$

对于 $A(x),B(x)$ 我们任取 $n+m+1$ 个点, $(x_0,A(x_0)B(x_0))...(x_{n+m+1},A(x_{n+m+1})B(x_{n+m+1}))$ 这个东西我们可以 $O(n+m)$ 求出

即我们可以在 $O(n)$ 时间内求出用来表示 $C(x)$ 的点

有了这些点,我们还需要快速的进行系数表示法与点表示法之间的相互转化这个操作,这就是FFT的核心

对于多项式 $A(x)=a_0+a_1x+...+a_{n-1}x^{n-1}$,假设 $n=2^k$,我们将它的奇数项与偶数项分类

令 $A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1},A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1}$

则 $A(x)=A_1(x^2)+xA_2(x^2)$

考虑取复数的 $n$ 次单位根上的点,

对于 $k\in[0,\frac{n}{2}-1],A(x)=A_1(w_n^{2k})+w_n^kA_2(w_n^{2k})=A_1(w_{\frac{n}{2}}^k)+w_{\frac{n}{2}}^kA_2(w_{\frac{n}{2}}^k)$

对于 $k\in[\frac{n}{2},n-1]$,直接对上面的式子加上一个 $\frac{n}{2}$, $A(x)=A_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}A_2(w_n^{2k+n})=A_1(w_{\frac{n}{2}}^k)-w_{\frac{n}{2}}^kA_2(w_{\frac{n}{2}}^k)$

那么对于每一个区间 $A(x)$,我们都可以将它分成两个区间来求

一共只会有 $log(n)$ 层,每次的长度为 $n$,所以我们可以在时间复杂度为 $O(nlog(n))$ 内求出这些点

到此,这就是傅里叶变换的正变换

接下来考虑怎样将点表示法转换为系数表示法

同样令 $A(x)=a_0+a_1x+...+a_{n-1}x^{n-1}$,$(w_n^k,A(w_n^k)=y_k)$ 表示我们取得 $n$ 个点

令 $c_k=\sum\limits_{i=1}^{n-1}y_i(w_n^{-k})^i$,则有 $c_k=na_k$

证明

$$c_k=\sum\limits_{i=0}^{n-1}(w_n^{-k})^i\sum\limits_{j=0}^{n-1}a_j(w_n^i)^j\\=\sum\limits_{i=0}^{n-1}(w_n^{-k})^i\sum\limits_{j=0}^{n-1}a_j(w_n^j)^i\\=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j(w_n^{j-k})^i\\=\sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}(w_n^{j-k})^i$$

令 $f(x)=1+x+x^2+...+x^{n-1}$

则 $xf(x)=x+x^2+...+x^n\Rightarrow (1-x)f(x)=1-x^n\Rightarrow (1-w_n^k)f(w_n^k)=1-(w_n^k)^n=0$

1.当 $k\ne 0$ 时,$f(w_n^k)=0$

2.当 $k=0$ 时,$f(w_n^k)=f(1)=n$

所以 $c_k=na_k$

证毕

也就是说,我们只需要快速求出 $c_k$ 就行了

不妨设 $B(x)=y_0+y_1x+...+y_{n-1}x^{n-1}$

那么 $c_k=B(w_n^{-k})$,这不就和上面求 $A(x)$ 是一样的吗,求出 $c_k$, $a_k=\frac{c_k}{n}$

到此就是傅里叶变换的全部内容了

###迭代实现FFT

根据前人的经验总结,递归写FFT常数会很大,所以一般用迭代来写

如果要用迭代来写的话,那么我们就需要处理出最后一层的顺序,举个例子

如图

可以观察出最后一层的下标用二进制表示出来与第一层的下标似乎有某种关系

可以证明

对于一个二进制表示为 $1001...100?$, 若最后一位为 $1$,则会被划分到右边的区间,即高位为 $1$,反之划分到左边区间,即高位为 $0$,

我们可以用递推求出每一个位置在底层的数

即用我们要求的数先右移一位,用之前求得的该位置的数或上最高位

那么这样,我们就可以最后一层往上迭代,就可以用迭代来实现傅里叶变换了

P3803 【模板】多项式乘法（FFT）

结合代码注释进行理解

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const double PI = acos(-1);
const int N = 3e5 + 5;

struct Complex
{
    double x, y;
    Complex operator + (const Complex t) const
    {return {x + t.x, y + t.y};}
    Complex operator - (const Complex t) const
    {return {x - t.x, y - t.y};}
    Complex operator * (const Complex t) const
    {return {x * t.x - y * t.y, x * t.y + y * t.x};}
} a[N], b[N];//结构体定义复数运算
int rev[N], bit, tot;

void FFT(Complex a[], int inv)
{
    for (int i = 0; i < tot; i++)
        if (i < rev[i])
            swap(a[i], a[rev[i]]);//求出低层所有的数,加判断语句防止多次交换
    for (int mid = 1; mid < tot; mid <<= 1)//枚举区间长度
    {
        Complex w1 = Complex({cos(PI / mid), inv * sin(PI / mid)});
        for (int i = 0; i < tot; i += mid * 2)//枚举每一个区间
        {
            Complex wk = Complex({1, 0});
            for (int j = 0; j < mid; j++, wk = wk * w1)//往上迭代
            {
                Complex x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}

int main()
{
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; i++)
        scanf("%lf", &a[i].x);
    for (int i = 0; i <= m; i++)
        scanf("%lf", &b[i].x);//一开始读入只有实部
    while ((1 << bit) < n + m + 1) bit++;//找到大于n+m+1的二的整次幂,方便迭代
    tot = 1 << bit;
    for (int i = 0; i < tot; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));//预处理出所有的在最后一层的位置
    FFT(a, 1), FFT(b, 1);
    for (int i = 0; i < tot; i++)
        a[i] = a[i] * b[i];//点表示法算卷积
    FFT(a, -1);//逆变换
    for (int i = 0; i <= n + m; i++)
        printf("%d ", (int)(a[i].x / tot + 0.5));//输出,加0.5为了保证精度
    return 0;
}