多项式算法——快速傅里叶变换FFT

多项式算法——快速傅里叶变换FFT

两个多项式函数相加的时间复杂度为$O(n)$，而相乘的时间复杂度为$O(n^2)$，即所有系数的笛卡尔积。傅里叶变换DFT将优化多项式相乘（以下简称卷积），而快速傅里叶变换则优化至$O(n\log n)$。

前置知识：

1. 微积分（欧拉公式）
2. 复变函数（单位元根）
3. 线性代数（矩阵方程）

多项式

本章定义的多项式函数为定义在实数域上，关于变量$x$形如：

$$A(x)=\sum _{j=0}^{n-1}{a}_j{x}^j$$

的$n$次多项式，$a_j$为多项式的系数，如果多项式非零最高次数项的系数为$a_k$，那么我们说$A(x)$的次数为$k$，记作$degree(A)=k$。

任何一个严格大于多项式$A(x)$的次数的整数都是多项式$A(x)$的一个次数界。

多项式加法

两个次数界为$n$的多项式$A(x)$和$B(x)$的多项式相加，仍是一个次数界为$n$的多项式$C(x)$，点值满足$C(x)=A(x)+B(x)$，且系数满足：

$$C(x)=\sum _{j=0}^{n-1}(a_j+b_j)x^j$$

多项式乘法

两个次数界为$n$的多项式$A(x)$和$B(x)$的多项式相乘，其结果是一个次数界为$2n-1$的多项式$C(x)$。点值满足$C(x)=A(x)B(x)$，且系数满足：

$$C(x)=\sum _{j=0}^{2n-1}\sum _{i=0}^j{a}_i{b}_{j-i}{x}^j$$

也称卷积。

多项式表示法

系数表示

对于一个次数为$n$的多项式$A(x)={\sum }_{j=0}^{n-1}{a}_j{x}^j$而言，其系数表示法是一组由系数组成的$n$维列向量$\vec{a}=(a_0,a_1,\dots ,a_{n-1})$。

对于给定的一个系数表示法的多项式，我们可以通过霍纳法则（秦九韶算法）在$O(n)$时间内计算出点值$A(x_0)$，即：

$$A(x_0)=a_0+x_0(a_1+x_0(a_2+\dots +x_0(a_{n-1}))\dots ))$$

同时，我们对两个系数表示法的多项式只需要$O(n)$的时间内就能完成，只需要做一次向量加法即可。

现在考虑系数表示法的多项式相乘，需要在$O(n^2)$的时间内完成，因为需要做一次向量卷积（笛卡尔直积），记作$\vec{c}=\vec{a}\ast \vec{b}$。

点值表示

一个次数界为$n$的多项式的点值表示就是由一个$n$个点值对所组成的集合：

$$\{(x_0,y_0),\dots ,(x_{n-1},y_{n-1})\}$$

$y_k=A(x_k)$，使得所有取值点的$x_k$均不相同。

一个多项式的取值点的不同，表达也不同，故一个多项式的点值表示并不唯一。我们统一取值点为一个固定的点集$x_0,\dots ,x_{n-1}$，称为基底。

求值计算的逆（从点值表示法计算系数表示法的过程）叫做插值。

定理：只有点值集合的大小和多项式的次数界相同，其插值才是唯一的。

证明：

点值和插值可以由如下矩阵表示：

$$\left[\begin{array}{cccc}1& x_0& \cdots & x_0^{n-1}\\ 1& x_1& \cdots & x_1^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{n-1}& \cdots & x_{n-1}^{n-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n-1\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n-1\end{array}\right]$$

左边的矩阵为范德蒙德矩阵，其行列式的值为${\prod }_{0\le j<k\le n-1}(x_k-x_j)$。

因为不存在两个相同的$x_k$，因此该矩阵可逆，存在逆矩阵，该矩阵方程有唯一解。

该方法算作一个插值方法，通过LU等分解方式，加速矩阵乘积运算，可在$O(n^3)$时间内进行插值。

作为了解，一种更快的插值方式是使用拉格朗日插值公式，其可以在$O(n^2)$内完成插值。

因此，求值和插值是定义完备的互逆运算。

我们可以看到，通过点值计算具有相同基底多项式乘法的时间复杂度为$O(n)$，即基底不变$y_k$对应相乘。同理，多项式加法的时间复杂度也为$O(n)$。

计算一个点值的时间复杂度是$O(n)$的，计算$n$个不同的点的时间复杂度是$O(n^2)$的，我们想法巧妙的构造一组基底，使得优化到$O(n\log n)$。

我们通过选取$2n$个单位复数根（记作${\omega }_{2n}$）作为基底，可以将求值和插值这两个过程优化到$O(n\log n)$。

其通过FFT计算多项式乘法定义了四个过程：

1. 加倍次数界：因为结果的多项式的界是$2n$的，因此我们需要取$2n$个点。
2. 求值：利用FFT算法，对两个多项式进行求值。
3. 逐点相乘：对应点相乘。
4. 插值：利用逆FFT算法，对多项式进行插值。

1和3步的时间为$O(n)$，24步的时间为$O(n\log n)$。

DFT与FFT

单位复数根

n次单位复数根是满足${\omega }^n=1$的复数$\omega$，这样的不同的复数恰好有$n$个，对于$k=0,\dots ,n-1$，$\omega =e^{\frac{2k\pi }{n}i}$。为了解释这个表达式，我们利用复数指数的形式解释：

$$e^{iu}=\cos (u)+i\sin (u)$$

故${\omega }^n=e^{2k\pi i}=\cos (2k\pi )+i\sin (2k\pi )=1$，其中指数$\frac{2k\pi }{n}$，恰好把一个单位圆周分成$n$分，正好对应着$n$个不同的解，当$k=1$时，${\omega }_n=e^{\frac{2\pi i}{n}}$，称为主$n$次单位复根，其他解都是主$n$次单位复根的幂次，故下文表示$n$个不同的根时，用${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$表示。

一般的，n次单位负数根在乘法意义上形成一个群，因为三角函数的周期性，n次单位复数根也存在周期性。

关于单位复数根的引理和证明

消去引理

对于任何的整数$n\ge 0,k\ge 0,d>0$，有：

$${\omega }_{dn}^{dk}=e^{\frac{2dk\pi }{dn}i}=e^{\frac{2k\pi }{n}i}={\omega }_n^k$$

推论：

$${\omega }_n^{\frac{n}{2}}={\omega }_2=-1$$

折半引理

如果$n>0$是偶数，那么$n$个$n$次单位复数根的平方的集合，等同于$n/2$个$n/2$次单位复数根的的集合。

根据消去引理$({\omega }_n^k{)}^2={\omega }_{\frac{n}{2}}^k$，如果对每个根都进行平方，那么每个不同的数正好出现两次，因为：

$$({\omega }_n^{k+n/2}{)}^2={\omega }_n^{2k+n}={\omega }_n^{2k}{\omega }_n^n={\omega }_n^{2k}=({\omega }_n^k{)}^2={\omega }_{n/2}^k$$

因此，$({\omega }_n^{k+n/2}{)}^2$和$({\omega }_n^k{)}^2$和${\omega }_{n/2}^k$相同。

注意，折半引理对于分治FFT是十分有帮助的，因为他保证了子问题的父问题的规模的一半。

求和引理

对于任意的整数$n>0$，和不能被$n$整除的非负整数$k$，存在：

$$\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=0$$

证明使用等比数列求和公式即可，这里不再证明。

DFT

我们希望在${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$处计算$A({\omega }_n^k)$的值，即在$n$次单位复数根处进行求值。

如果多项式$A(x)={\sum }_{j=0}^{n-1}{a}_j{x}^j$，那么$y_k=A({\omega }_n^k)={\sum }_{j=0}^{n-1}{a}_j{\omega }_n^{kj}$。

那么向量$\vec{y}=(y_0,\dots ,y_{n-1})$就是系数向量$\vec{a}=(a_0,\dots ,a_{n-1})$的离散傅里叶变换（DFT）。记作$\vec{y}={\text{DFT}}_n(\vec{a})$。

FFT

这里介绍一种快速傅里叶变换（FFT）的算法，利用单位复数根的性质，我们就可在$O(n\log n)$的时间内计算出${\text{DFT}}_n(\vec{a})$，通篇假设$n$都是$2$的整数次幂。

FFT利用了分治策略，分别采用偶数下标和奇数下标的系数构造两个新的多项式：

$$\begin{gathered} A_0(x)=a_0+a_{2}x+a_4{x}^2+\dots +a_{n-2}{x}^{n/2-1} \\ {A}_1(x)=a_1+a_{3}x+a_5{x}^2+\dots +a_{n-1}{x}^{n/2-1} \end{gathered}$$

其中$A_0(x)$包含偶数下标，$A_1(x)$包含奇数下标，于是有：

$$A(x)=A_0(x^2)+x{A}_1(x^2)$$

所以，求$A$在${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$处的值就转换为求界为：

1. $n/2$的两个多项式$A_0$和$A_1$在点$({\omega }_n^0{)}^2,({\omega }_n^1{)}^2,\dots ,({\omega }_n^{n-1}{)}^2$处的取值
2. 根据$A(x)=A_0(x^2)+x{A}_1(x^2)$合并结果。

根据折半引理，如果$n>0$是偶数，那么$n$个$n$次单位复数根的平方的集合，等同于$n/2$个$n/2$次单位复数根的的集合。这样，问题就变成了两个界为$n/2$，单位复根数量为$n/2$的FFT问题，此时问题变成了两个子问题，但问题规模减半。

算法：

typedef complex<double> dcpx;

void recur_fft(dcpx A[], dcpx B[], int n, int inv)
{
    if (n == 1)
        return;
    int hf = n >> 1;
    for (int i = 0; i < hf; i++)
    {
        B[i] = A[i << 1];
        B[i + hf] = A[(i << 1) | 1];
    }
    for (int i = 0; i < n; i++)
        A[i] = B[i];

    recur_fft(A, B, hf, inv);
    recur_fft(A + hf, B + hf, hf, inv);

    dcpx omega_n(cos(2 * M_PI / n), inv * sin(2 * M_PI / n));
    dcpx omega(1, 0);
    for (int i = 0; i < hf; i++)
    {
        B[i] = A[i] + omega * A[i + hf];
        B[i + hf] = A[i] - omega * A[i + hf];
        omega *= omega_n;
    }
    for (int i = 0; i < n; i++)
        A[i] = B[i];
}

前几行按照奇偶性分类系数，然后递归，主要在最后循环合并答案上面，该代码$A$代表输入的系数向量，$B$代表辅助数组，运行的结果将变换的结果放在$A$向量里，参数$inv$稍后再说他的作用，在这里我们传入$1$。

循环的最后一行，omega *= omega_n;能保证我们能遍历$w_n^0,w_n^1,w_n^{n/2-1}$。

考虑$k<\frac{n}{2}$，此时表示成$k=i$，那么$y_k=A_0({\omega }_n^{2k})+{\omega }_n^k{A}_1({\omega }_n^{2k})=A_0({\omega }_{n/2}^k)+{\omega }_n^k{A}_1({\omega }_{n/2}^k)=y_0[i]+{\omega }_n^i{y}_1[i]$。

考虑$k\ge \frac{n}{2}$，此时表示成$k=i+n/2$，那么$y_k=A_0({\omega }_{n/2}^{i+n/2})+{\omega }_n^{i+n/2}{A}_1({\omega }_{n/2}^{i+n/2})=A_0({\omega }_{n/2}^i)-{\omega }_n^i{A}_1({\omega }_{n/2}^i)=y_0[i]-{\omega }_n^i{y}_1[i]$。

前面的系数$\pm {\omega }_n^i$，也称旋转因子。我们知道，乘以一个虚数相当于进行旋转。

根据主定理，其时间复杂度为$O(n\log n)$。

逆FFT

考虑复矩阵方程：

$$\left[\begin{array}{cccc}1& {\omega }_0& \cdots & {\omega }_0^{n-1}\\ 1& {\omega }_1& \cdots & {\omega }_1^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_{n-1}& \cdots & {\omega }_{n-1}^{n-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n-1\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n-1\end{array}\right]$$

记最左边的矩阵为$V_n$，我们的目标就是找到$V_n$的逆矩阵$V_n^{-1}$，进而求解复矩阵方程。

定理：对于$j,k$，$V_n[j,k]=\frac{{\omega }_n^{-jk}}{n}$。

证明：考虑证明$V_n{V}_n^{-1}=I_n$即可。需要使用上式求和引理。

给定逆矩阵$V_n^{-1}$，可以推导出$a_j$：

$$a_j=\frac{1}{n}\sum _{k=0}^{n-1}{y}_k{\omega }_n^{-kj}$$

我们和式子$y_k=A({\omega }_n^k)={\sum }_{j=0}^{n-1}{a}_j{\omega }_n^{kj}$进行对比，发现逆FFT过程只多了一个负号和除以$n$，我们可以稍微修改一下FFT就可以完成逆FFT，即只需要把$inv$传入$-1$，并在最终除以$n$即可。

离散傅里叶变换卷积定理：对于任意两个长度为$n$的向量$a$和$b$，其中$n$是$2$的$n$次幂，则

$$a\ast b={\text{DFT}}_{2n}^{-1}({\text{DFT}}_{2n}(a)\cdot {\text{DFT}}_{2n}(b))$$

其中向量$a$和$b$通过系数补零的方式扩大值$2n$即可。

代码实现

P3803

typedef complex<double> dcpx;

void recur_fft(dcpx A[], dcpx B[], int n, int inv)
{
    if (n == 1)
        return;
    int hf = n >> 1;
    for (int i = 0; i < hf; i++)
    {
        B[i] = A[i << 1];
        B[i + hf] = A[(i << 1) | 1];
    }
    for (int i = 0; i < n; i++)
        A[i] = B[i];

    recur_fft(A, B, hf, inv);
    recur_fft(A + hf, B + hf, hf, inv);

    dcpx omega_n(cos(2 * M_PI / n), inv * sin(2 * M_PI / n));
    dcpx omega(1, 0);
    for (int i = 0; i < hf; i++)
    {
        B[i] = A[i] + omega * A[i + hf];
        B[i + hf] = A[i] - omega * A[i + hf];
        omega *= omega_n;
    }
    for (int i = 0; i < n; i++)
        A[i] = B[i];
}
const int MAXT = 1 << 21;
dcpx A[MAXT];
dcpx B[MAXT];
dcpx T[MAXT];
int main()
{
    FR;
    int n, m;
    scanf("%d %d", &n, &m);

    int mx = max(n, m);

    for (int i = 0; i < n + 1; i++)
    {
        int coff;
        scanf("%d", &coff);
        A[i].real(coff);
    }
    for (int i = 0; i < m + 1; i++)
    {
        int coff;
        scanf("%d", &coff);
        B[i].real(coff);
    }

    recur_fft(A, T, MAXT, 1);
    recur_fft(B, T, MAXT, 1);

    for (int i = 0; i < MAXT; i++)
    {
        A[i] *= B[i];
    }

    recur_fft(A, T, MAXT, -1);

    for (int i = 0; i < n + m + 1; i++)
    {
        ll ans = A[i].real() / MAXT + 0.5;
        printf("%lld ", ans);
    }
    return 0;
}

改进FFT

迭代版本的FFT

考虑我们的递归版的递归过程，我们是将其进行系数向量的拆解，然后逐层合并。如果我们可以在一开始就将系数向量按位置拆解，那么我们就可以通过循环的方法自底向上的方式进行合并答案。

考虑$n=8$的情况，重新排列系数下标依次为$a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7$。我们看看$0,4,2,6,1,5,3,7$和$0,1,2,3,4,5,6,7$的关系，我们发现两个数列的二进制每一项都是对应项的转置。这种方法叫做逆位置对换。

我们可以通过递推的方式来实现：

for (int i = 1; i < MAXT; i++)
{
    rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (E - 1));
}

考虑逆位置对换算法，我们遍历每一个位置$i$，如果$i<rev[i]$，说明$i$应该交换，且一次交换正好能归位两个元素，因为$rev[rev[i]]=i$。

for (int i = 0; i < n; i++)
    if (i < rev[i])
        swap(A[i], A[rev[i]]);

完整代码：

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

#define FR freopen("in.txt", "r", stdin)

typedef long long ll;

typedef complex<double> dcpx;

const int E = 21;
const int MAXT = 1 << E;
const double PI = acos(-1);

int rev[MAXT];

void iter_fft(dcpx A[], int LOG, int inv)
{
    int n = 1 << LOG;
    for (int i = 0; i < n; i++)
        if (i < rev[i])
            swap(A[i], A[rev[i]]);

    for (int e = 1; e <= LOG; e++)
    {
        int m = 1 << e;
        for (int i = 0; i < n; i += m)
        {
            int hf = m / 2;
            dcpx omega(1, 0);
            dcpx omega_n(cos(2 * PI / m), inv * sin(2 * PI / m));
            for (int j = 0; j < hf; j++)
            {
                dcpx x = A[i + j];
                dcpx y = A[i + j + hf] * omega;
                A[i + j] = x + y;
                A[i + j + hf] = x - y;
                omega *= omega_n;
            }
        }
    }
}

dcpx A[MAXT];
dcpx B[MAXT];

int main()
{
    for (int i = 1; i < MAXT; i++)
    {
        rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (E - 1));
    }
    int n, m;

    cin >> n >> m;

    for (int i = 0; i < n + 1; i++)
    {
        double a;
        cin >> a;
        A[i].real(a);
    }
    for (int i = 0; i < m + 1; i++)
    {
        double a;
        cin >> a;
        B[i].real(a);
    }

    iter_fft(A, E, 1);
    iter_fft(B, E, 1);

    for (int i = 0; i < MAXT; i++)
    {
        A[i] *= B[i];
    }

    iter_fft(A, E, -1);

    for (int i = 0; i < m + n + 1; i++)
    {
        cout << ll(A[i].real() / MAXT + 0.5) << " ";
    }
    return 0;
}

例题

P1919

高精度乘法，考虑一个数字的按权展开的形式：

$$a=\sum _{i=0}^{\dots }{a}_i\times 10^i$$

这是一个多项式，那么两个数的乘积就是两个多项式相乘，此时就可以使用FFT卷起来了。

最后在考虑进位问题即可。

const int E = 21;
const int MAXT = 1 << E;
int rev[MAXT];

void iter_fft(dcpx A[], int inv)
{
    for (int i = 0; i < MAXT; i++)
    {
        if (i < rev[i])
            swap(A[i], A[rev[i]]);
    }

    for (int e = 1; e <= E; e++)
    {
        int m = 1 << e;
        int hf = m >> 1;
        for (int i = 0; i < MAXT; i += m)
        {
            dcpx omega(1, 0);
            dcpx omega_n(cos(2 * M_PI / m), inv * sin(2 * M_PI / m));
            for (int k = 0; k < hf; k++)
            {
                dcpx x = A[i + k];
                dcpx y = A[i + k + hf] * omega;
                A[i + k] = x + y;
                A[i + k + hf] = x - y;
                omega *= omega_n;
            }
        }
    }
}

dcpx A[MAXT];
dcpx B[MAXT];

int main()
{
    FR;
    for (int i = 1; i < MAXT; i++)
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (E - 1));
    }
    string a, b;
    cin >> a >> b;
    for (int i = a.size() - 1; i >= 0; i--)
    {
        A[a.size() - i - 1] = a[i] - '0';
    }
    for (int i = b.size() - 1; i >= 0; i--)
    {
        B[b.size() - i - 1] = b[i] - '0';
    }

    iter_fft(A, 1);
    iter_fft(B, 1);

    for (int i = 0; i < MAXT; i++)
    {
        A[i] *= B[i];
    }

    iter_fft(A, -1);
    string ans;
    ll carry = 0;
    for (int i = 0; i < MAXT; i++)
    {
        ll pp = ll(A[i].real() / MAXT + 0.5) + carry;
        carry = pp / 10;
        ans.push_back((pp % 10) + '0');
    }
    int i = MAXT - 1;
    while (i >= 0 && ans[i] == '0')
        i--;
    string Ans;
    for (; i >= 0; i--)
    {
        Ans.push_back(ans[i]);
    }
    cout << Ans;
    return 0;
}

总结

本文介绍了多项式的基础和FFT算法的实现，是多项式算法系列的基础。