浅谈算法——从多项式乘法到FFT

多项式乘法

我们知道，多项式可以表示成：

A=∑i=0naixi$$A=\sum _{i=0}^n{a}_i{x}^i$$

的形式。
对于两个多项式A(x)$A(x)$和B(x)$B(x)$，我们可以计算乘积A⋅B$A\cdot B$：

A⋅B=∑i=0sizeA∑j=0sizeBaibjxi+j$$A\cdot B=\sum _{i=0}^{siz{e}_A}\sum _{j=0}^{siz{e}_B}{a}_i{b}_j{x}^{i+j}$$

但是，这样算是O(sizeA⋅sizeB)$O(siz{e}_A\cdot siz{e}_B)$的，太慢了，怎么办？
我们需要换一条思路。

首先，我们得知道一个东西：多项式的点值表示法。
我们把上面的称为多项式的系数表示法，而点值表示法就是：
若A$A$多项式的次数为n$n$，则任取n$n$个不相同的x0,x1,⋯,xn$x_0,x_1,\cdots ,x_n$，求出A$A$多项式的A(x0),A(x1),⋯,A(xn)$A(x_0),A(x_1),\cdots ,A(x_n)$。记为：

<(x0,A(x0)),(x1,A(x1)),⋯,(xn,A(xn))>$$<(x_0,A(x_0)),(x_1,A(x_1)),\cdots ,(x_n,A(x_n))>$$

显然，一个有n+1$n+1$个点的点对唯一表示一个n$n$次多项式。

对于一个点值表示法下多项式

<(x0,A(x0)),(x1,A(x1)),⋯,(xn,A(xn))>$$<(x_0,A(x_0)),(x_1,A(x_1)),\cdots ,(x_n,A(x_n))>$$

和

<(x0,B(x0)),(x1,B(x1)),⋯,(xn,B(xn))>$$<(x_0,B(x_0)),(x_1,B(x_1)),\cdots ,(x_n,B(x_n))>$$

它们的乘积是

<(x0,A(x0)⋅B(x0)),(x1,A(x1)⋅B(x1)),⋯,(xn,A(xn)⋅B(xn))>$$<(x_0,A(x_0)\cdot B(x_0)),(x_1,A(x_1)\cdot B(x_1)),\cdots ,(x_n,A(x_n)\cdot B(x_n))>$$

可以看出点值表示法的多项式相乘是O(n)$O(n)$的。

等等，我们好像找到了一个突破口！
为啥不把原来的系数表示法下的多项式转化成点值表示法呢？
仔细想一想：系数表示法与点值表示法互相转换，这个步骤好像是O(n2)$O(n^2)$的。
而FFT（快速傅里叶变换）就是为了优化这个O(n2)$O(n^2)$。

PS：对于O(n2)$O(n^2)$的点值表示法转化成系数表示法可以看百度百科中关于插值法的介绍。

FFT（快速傅里叶变换）

如果未特别说明，那么下面的多项式次数将是2k−1$2^k-1$的形式。
如果不是关键部分的公式或定理，不提供证明，自己出门右转百度。

首先介绍两个概念：
DFT（离散傅里叶变换）是将多项式由系数表示法转化成点值表示法；
IDFT（离散傅里叶逆变换）是将多项式由点值表示法转化成系数表示法；
而FFT就是上述两种变换的优化。

DFT部分

前置技能：
下面的内容将会提到复数，不会的可以参考百度百科中关于复数的介绍；

为了介绍FFT中的DFT部分，首先要介绍的是一个概念：单位根。
单位根：若有

zn=1$$z^n=1$$

此时将z$z$称为n$n$次单位根。
若有z∈R$z\in \mathbb{R}$，显然，z$z$可以等于1$1$，如果n$n$是偶数，则z$z$还可以等于−1$-1$。
我们把范围扩大到z∈C$z\in \mathbb{C}$，那么，我们可以得到n$n$个复数，它们将复平面上的单位圆等分成n$n$份。

为了表示n$n$次单位根，我们引入一个公式。
欧拉公式：

ein=cosn+isinn$$e^{in}=\cos n+i\sin n$$

如果我们令：

ωn=e2πi/n$${\omega }_n=e^{2\pi i/n}$$

那么，n$n$次单位根就可以表示成ω0n,ω1n,⋯,ωn−1n${\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}$，它们的n$n$次方显然都是1$1$。

下面是关于ωn${\omega }_n$的两条性质：（都是在n$n$为偶数的情况下）

ωn/2n=e(2πi/n)⋅(n/2)=eπi=cosπ+isinπ=−1(1)$$\begin{array}{}\text{(1)}& {\omega }_n^{n/2}=e^{(2\pi i/n)\cdot (n/2)}=e^{\pi i}=cos\pi +i\sin \pi =-1\end{array}$$

ω2n=e2⋅2πi/n=ωn/2(2)$$\begin{array}{}\text{(2)}& {\omega }_n^2=e^{2\cdot 2\pi i/n}={\omega }_{n/2}\end{array}$$

下面，我们进入正题：DFT的求法。
在这里，我们令多项式次数为n−1$n-1$，那么我们可以用点值表示成

<(ω0n,A(ω0n),(ω1n,A(ω1n)),⋯,(ωn−1n,A(ωn−1n))>$$<({\omega }_n^0,A({\omega }_n^0),({\omega }_n^1,A({\omega }_n^1)),\cdots ,({\omega }_n^{n-1},A({\omega }_n^{n-1}))>$$

额……这时间复杂度好像并没有减少……
别急，我们来看A(ωkn)$A({\omega }_n^k)$能够表示成什么。

A(ωkn)====∑i=0n−1aiωkin∑i=0n/2−1a2iω2kin+∑i=0n/2−1a2i+1ω2ki+kn∑i=0n/2−1a2iωkin/2+∑i=0n/2−1a2i+1ωkin/2ωkn∑i=0n/2−1a2iωkin/2+ωn∑i=0n/2−1a2i+1ωkin/2(3)(4)(5)(6)$$\begin{array}{}\text{(3)}& A({\omega }_n^k)=& \sum _{i=0}^{n-1}{a}_i{\omega }_n^{ki}\text{(4)}& =& \sum _{i=0}^{n/2-1}{a}_{2i}{\omega }_n^{2ki}+\sum _{i=0}^{n/2-1}{a}_{2i+1}{\omega }_n^{2ki+k}\text{(5)}& =& \sum _{i=0}^{n/2-1}{a}_{2i}{\omega }_{n/2}^{ki}+\sum _{i=0}^{n/2-1}{a}_{2i+1}{\omega }_{n/2}^{ki}{\omega }_n^k\text{(6)}& =& \sum _{i=0}^{n/2-1}{a}_{2i}{\omega }_{n/2}^{ki}+{\omega }_n\sum _{i=0}^{n/2-1}{a}_{2i+1}{\omega }_{n/2}^{ki}\end{array}$$

我们来分别看一看这神奇的步骤。
(3)$(3)$步骤就是将ωkn${\omega }_n^k$带入原来的A$A$多项式。
(4)$(4)$步骤就是将原多项式拆成两个部分，按奇偶分类。
(5)$(5)$步骤用到了上面提到的性质(2)$(2)$。
(6)$(6)$步骤就是上面式子的后半部分提出公因数。

有了这个等式，我们就可以分治+递归解决DFT了。
算法步骤：

• 对当前的多项式（一个数组）系数进行奇偶分类；
• 递归算出偶数部分的数组的anse$an{s}_e$和奇数部分的数组的anso$an{s}_o$；
• 这个多项式的ans=anse+ωnanso$ans=an{s}_e+{\omega }_{n}an{s}_o$

但是这个的常数好像很大啊？能不能减少一点呢？
上面的性质(1)$(1)$给了我们提示：

ωn/2+kn=ωn/2n⋅ωkn=−ωkn$${\omega }_n^{n/2+k}={\omega }_n^{n/2}\cdot {\omega }_n^k=-{\omega }_n^k$$

在算k<n2$k<\frac{n}{2}$时，可以顺便把k≥n2$k\ge \frac{n}{2}$的情况也算出来。

常数减小了一半！但是还是很大啊！
递归版的程序一般比非递归版慢，为啥不用非递归版呢？

算法核心就是奇偶分类，分来分去最后分到了哪里？我们来研究研究。
显然，一个序列原来是0,1,2,3,4,5,6,7$0,1,2,3,4,5,6,7$，最终变成0,4,2,6,1,5,3,7$0,4,2,6,1,5,3,7$。
把它们的二进制列出来：

000,001,010,011,100,101,110,111000,100,010,110,001,101,011,111$$\begin{gathered} 000,001,010,011,100,101,110,111 \\ 000,100,010,110,001,101,011,111 \end{gathered}$$

其中，上面是位置，下面是这个位置对应的数。

把上面的数翻转，好像就是下面的数！
没错，只需要计算一下每个数的二进制翻转后的结果，就能得到一个数最终对应的位置，也就能实现非递归版了。

代码：

int dft_fast(complex* ar,int len)
{
  for(register int i=0; i<len; ++i)
    {
      if(rev[i]<i)
        {
          std::swap(ar[rev[i]],ar[i]);//交换一个位置和它的翻转后位置
        }
    }
  for(register int i=2; i<=len; i<<=1)//i代表当前序列的长度
    {
      complex wn(cos(2*pi/i),sin(2*pi/i));//omega_n
      for(register int j=0; j<len; j+=i)//j代表序列的起始位置
        {
          complex w(1,0);//下面代表omega_n^k
          for(register int k=0; k<(i>>1); ++k)//枚举i次单位根的每一种取值
            {
              complex x=ar[j+k],y=w*ar[j+k+(i>>1)];
              ar[j+k]=x+y;//合并操作，将两边合并成一个点值表示法下的多项式
              ar[j+k+(i>>1)]=x-y;
              w=w*wn;
            }
        }
    }
  return 0;
}

IDFT部分

回顾上面的DFT部分，仔细思考一下，它本质就是在求：

⎧⎩⎨⎪⎪⎪⎪⎪⎪a0(ω0n)0+a1(ω0n)1+⋯+an−1(ω0n)n−1=A(ω0n)a0(ω1n)0+a1(ω1n)1+⋯+an−1(ω1n)n−1=A(ω1n)⋯a0(ωn−1n)0+a1(ωn−1n)1+⋯+an−1(ωn−1n)n−1=A(ωn−1n)$$\{\begin{array}{l}{a}_0({\omega }_n^0{)}^0+a_1({\omega }_n^0{)}^1+\cdots +a_{n-1}({\omega }_n^0{)}^{n-1}=A({\omega }_n^0)\\ a_0({\omega }_n^1{)}^0+a_1({\omega }_n^1{)}^1+\cdots +a_{n-1}({\omega }_n^1{)}^{n-1}=A({\omega }_n^1)\\ \cdots \\ a_0({\omega }_n^{n-1}{)}^0+a_1({\omega }_n^{n-1}{)}^1+\cdots +a_{n-1}({\omega }_n^{n-1}{)}^{n-1}=A({\omega }_n^{n-1})\end{array}$$

其中，给定了a0,a1,⋯,an−1$a_0,a_1,\cdots ,a_{n-1}$以及ω0n,ω1n,⋯,ωn−1n${\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}$，求A(ω0n),A(ω1n),⋯,A(ωn−1n)$A({\omega }_n^0),A({\omega }_n^1),\cdots ,A({\omega }_n^{n-1})$的值。

用矩阵表示如下：

⎡⎣⎢⎢⎢⎢⎢(ω0n)0(ω1n)0⋮(ωn−1n)0(ω0n)1(ω1n)1⋮(ωn−1n)1⋯⋯⋱⋯(ω0n)n−1(ω1n)n−1⋮(ωn−1n)n−1⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢A(ω0n)A(ω1n)⋮A(ωn−1n)⎤⎦⎥⎥⎥⎥⎥(7)$$\begin{array}{}\text{(7)}& \left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{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}A({\omega }_n^0)\\ A({\omega }_n^1)\\ ⋮\\ A({\omega }_n^{n-1})\end{array}\right]\end{array}$$

我们令：

V=⎡⎣⎢⎢⎢⎢⎢(ω0n)0(ω1n)0⋮(ωn−1n)0(ω0n)1(ω1n)1⋮(ωn−1n)1⋯⋯⋱⋯(ω0n)n−1(ω1n)n−1⋮(ωn−1n)n−1⎤⎦⎥⎥⎥⎥⎥$$V=\left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right]$$

那么IDFT的本质就是求V$V$矩阵的逆矩阵。

考虑下面这个矩阵：

D=⎡⎣⎢⎢⎢⎢⎢⎢(ω−0n)0(ω−1n)0⋮(ω−(n−1)n)0(ω−0n)1(ω−1n)1⋮(ω−(n−1)n)1⋯⋯⋱⋯(ω−0n)n−1(ω−1n)n−1⋮(ω−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥⎥$$D=\left[\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& \cdots & ({\omega }_n^{-0}{)}^{n-1}\\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& \cdots & ({\omega }_n^{-1}{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{-(n-1)}{)}^0& ({\omega }_n^{-(n-1)}{)}^1& \cdots & ({\omega }_n^{-(n-1)}{)}^{n-1}\end{array}\right]$$

那么我们令E=D⋅V$E=D\cdot V$，则：

Ei,j===∑k=0n−1Di,kVk,j∑k=0n−1(ω−in)k(ωkn)j∑k=0n−1ωk(j−i)n$$\begin{array}{rl}{E}_{i,j}=& \sum _{k=0}^{n-1}{D}_{i,k}{V}_{k,j}\\ =& \sum _{k=0}^{n-1}({\omega }_n^{-i}{)}^k({\omega }_n^k{)}^j\\ =& \sum _{k=0}^{n-1}{\omega }_n^{k(j-i)}\end{array}$$

显然：

Ei,j={0(i≠j)n(i=j)$$E_{i,j}=\{\begin{array}{l}0(i\ne j)\\ n(i=j)\end{array}$$

因此，

1nD⋅V=1nE=In$$\frac{1}{n}D\cdot V=\frac{1}{n}E=I_n$$

。
(7)$(7)$式两边同时左乘1nD$\frac{1}{n}D$，可得

⎡⎣⎢⎢⎢⎢a0a1⋮an−1⎤⎦⎥⎥⎥⎥=1n⎡⎣⎢⎢⎢⎢⎢⎢(ω−0n)0(ω−1n)0⋮(ω−(n−1)n)0(ω−0n)1(ω−1n)1⋮(ω−(n−1)n)1⋯⋯⋱⋯(ω−0n)n−1(ω−1n)n−1⋮(ω−(n−1)n)n−1⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢A(ω0n)A(ω1n)⋮A(ωn−1n)⎤⎦⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\frac{1}{n}\left[\begin{array}{cccc}({\omega }_n^{-0}{)}^0& ({\omega }_n^{-0}{)}^1& \cdots & ({\omega }_n^{-0}{)}^{n-1}\\ ({\omega }_n^{-1}{)}^0& ({\omega }_n^{-1}{)}^1& \cdots & ({\omega }_n^{-1}{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{-(n-1)}{)}^0& ({\omega }_n^{-(n-1)}{)}^1& \cdots & ({\omega }_n^{-(n-1)}{)}^{n-1}\end{array}\right]\left[\begin{array}{c}A({\omega }_n^0)\\ A({\omega }_n^1)\\ ⋮\\ A({\omega }_n^{n-1})\end{array}\right]$$

这就相当于把DFT中ωkn${\omega }_n^k$都换成ω−kn${\omega }_n^{-k}$。

FFT总代码

int fft(complex* ar,int len,int op)
{
  for(register int i=0; i<len; ++i)
    {
      if(rev[i]<i)
        {
          std::swap(ar[rev[i]],ar[i]);
        }
    }
  for(register int i=2; i<=len; i<<=1)
    {
      complex wn(cos(2*pi/i),sin(2*pi*op/i));//只有这里较DFT代码有变动
      for(register int j=0; j<len; j+=i)
        {
          complex w(1,0);
          for(register int k=0; k<(i>>1); ++k)
            {
              complex x=ar[j+k],y=w*ar[j+k+(i>>1)];
              ar[j+k]=x+y;
              ar[j+k+(i>>1)]=x-y;
              w=w*wn;
            }
        }
    }
  if(op==-1)
    {
      for(register int i=0; i<len; ++i)
        {
          ar[i].r/=len;
        }
    }
  return 0;
}

多项式乘法模板

#include <cstdio>
#include <cmath>
#include <algorithm>

const int maxn=100000;
const double pi=acos(-1);

struct complex
{
  double r,i;

  complex(double r_=0,double i_=0)
  {
    r=r_;
    i=i_;
  }

  complex operator +(const complex &other)
  {
    return complex(r+other.r,i+other.i);
  }

  complex operator -(const complex &other)
  {
    return complex(r-other.r,i-other.i);
  }

  complex operator *(const complex &other)
  {
    return complex(r*other.r-i*other.i,r*other.i+i*other.r);
  }
};

complex a[maxn<<2],b[maxn<<2],c[maxn<<2];
int rev[maxn<<2],n,m;

int fft(complex* ar,int len,int op)
{
  for(register int i=0; i<len; ++i)
    {
      if(rev[i]<i)
        {
          std::swap(ar[rev[i]],ar[i]);
        }
    }
  for(register int i=2; i<=len; i<<=1)
    {
      complex wn(cos(2*pi/i),sin(2*pi*op/i));
      for(register int j=0; j<len; j+=i)
        {
          complex w(1,0);
          for(register int k=0; k<(i>>1); ++k)
            {
              complex x=ar[j+k],y=w*ar[j+k+(i>>1)];
              ar[j+k]=x+y;
              ar[j+k+(i>>1)]=x-y;
              w=w*wn;
            }
        }
    }
  if(op==-1)
    {
      for(register int i=0; i<len; ++i)
        {
          ar[i].r/=len;
        }
    }
  return 0;
}

int main()
{
  scanf("%d%d",&n,&m);
  for(register int i=0; i<=n; ++i)
    {
      scanf("%lf",&a[i].r);
    }
  for(register int i=0; i<=m; ++i)
    {
      scanf("%lf",&b[i].r);
    }
  n=n+m;
  int l=0;
  m=1;
  while(m<=n)
    {
      ++l;
      m<<=1;
    }
  for(register int i=0; i<m; ++i)
    {
      rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    }
  fft(a,m,1);
  fft(b,m,1);
  for(register int i=0; i<m; ++i)
    {
      c[i]=a[i]*b[i];
    }
  fft(c,m,-1);
  for(register int i=0; i<n; ++i)
    {
      printf("%d ",(int)(c[i].r+0.5));
    }
  printf("%d\n",(int)(c[n].r+0.5));
  return 0;
}