avatar

目录
浅谈多项式与FFT

失智博主的学习日常之多项式与FFT
大部分摘自《算法导论》

多项式的表示

系数表达

  对于多项式 $A(x)=\sum\limits_{j=0}^{n-1}{a_j}{x^j}$ ,其系数表达就是由系数组成的向量 $a=(a_0,\ a_1,\ \cdots,\ a_{n-1})$ 。
  系数表达的多项式乘法 $C(x)=A(x)B(x)$ (注: $A$ 和 $B$ 在相同的 $n$ 位置取值),可以表达如下(也就是一项一项相乘,类似高精度):
$$C(x)=\sum\limits_{j=0}^{2n-2}{c_j x^{j}}$$

其中
$$c_j=\sum\limits_{k=0}^{j}{a_k b_{j-k}}$$

顺便一提,由上述公式得到的系数向量 $c$ 也称为系数向量 $a$ 和 $b$ 的卷积,记为 $c=a\otimes b$ 。

点值表达

  对于次数界(注:任何严格大于一个多项式次数的整数都是该多项式的次数界)为 $n$ 的多项式 $A(x)$ 的点值表达就是一个由 $n$ 个点值对所组成的集合$$\lbrace (x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\rbrace$$使得对 $k=0,\ 1,\ \cdots,\ n-1$ ,所有 $x_k$ 各不相同, $y_k=A(x_k)$ 。一个多项式可以有很多不同的点值表达。
  求值计算的逆(从一个多项式的点值表达式确定其系数表达式)成为插值。由插值多项式的唯一性(见下)可知,当插值多项式的次数界等于已知的点值对的数目,插值才是明确的。
  插值定理的唯一性:对于任意 $n$ 个点值对组成的集合 $\lbrace (x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\rbrace$ ,其中所有的 $x_k$ 都不同;那么存在唯一的次数界为 $n$ 的多项式 $A(x)$ ,满足 $y_k=A(x_k)$ ,$k=0,\ 1,\ \cdots,\ n-1$ 。
  相对系数表达,点值表达在多项式乘法方面是很方便的。如果 $C(x)=A(x)B(x)$ (注: $A$ 和 $B$ 在相同的 $n$ 位置取值),则对于任一点 $x_k$ ,$C(x_k)=A(x_k)B(x_k)$ 。但注意,若 $A$ 与 $B$ 的次数界为 $n$ ,此时 $C(x)$ 的次数界为 $2n$ 。则我们需要 $2n$ 个点对。
$$A:\lbrace (x_0,y_0),(x_1,y_1),\cdots,(x_{2n-1},y_{2n-1})\rbrace$$

$$B:\lbrace (x_0,y_0^{‘}),(x_1,y_1^{‘}),\cdots,(x_{2n-1},y_{2n-1}^{‘})\rbrace$$

$$C:\lbrace (x_0,y_0 y_0^{‘}),(x_1,y_1 y_1^{‘}),\cdots,(x_{2n-1},y_{2n-1} y_{2n-1}^{‘})\rbrace$$

  时间复杂度只有 $O(n)$ 。比采用系数表达法的 $O(n^2)$ 稳了很多。

系数形式表示的多项式快速乘法

  我们希望利用基于点值形式的表达式去加速基于系数形式表达的多项式乘法。那么关键就在于快速把一个多项式从系数形式转化为点值形式(求值),和从点值形式快速转化为系数形式(插值)。
  如果选取的点比较巧妙,则可以把两种表示法的转化降到 $O(nlog_2n)$ 级别。也就是接下来要说的。

DFT与FFT

单位复数根

单位复数根的定义

  如果随便选 $n$ 个点的话,暴力计算 $x_k^{0}, x_k^{1}, \cdots , x_k^{n-1}$ 的话需要 $O(n^2)$ 的时间。
  换种思路,我们带入一些特殊的 $x$ ,使得 $x$ 的若干次方等于1,我们就不用做次方运算了。我们需要的这种 $x$ 便是单位复数根
  我们称满足 $\omega ^{n}=1$ 的复数 $\omega$ 为 $n$ 次单位复数根。但是需要找的单位根这么多,仅实数的也就-1,1。
  我们需要一个非常优美的公式:欧拉公式

$$e^{ix}=cos x+isinx$$

当 $x=\pi$ 时,式子就成了这样:

$$e^{i \pi}+1=0$$

我们需要的是等于 $1$ 的,那就移项再平方一下:

$$e^{2\pi i}=1$$

然后我们就可以得到其中的一个 $n$ 次单位复数根了:

$$\omega _ {n}=e^{2\pi i/n}$$

这个单位根被称为主 $n$ 次单位根
  得到一个单位根后,我们可以对其进行幂运算,得到其他 $n-1$ 个 $n$ 次单位根。由于计算机能比较方便得计算三角函数所以结果我们一般用三角函数表示,证明时也会有 $e^n$ 形式的证明。
\begin{align}
\omega _ {n} ^ {k}&=(\omega _ {n})^{k} \\
&=e^{\frac{2k\pi i}{n}}\\
&=cos\frac{2k\pi}{n}+isin\frac{2k\pi}{n}\\
(k&=1,\ 2,\ \cdots,\ n-1)
\end{align}

  由三角形式的单位根表达式可以方便得看出 $n$ 次单位复数根一共有n个。对于大于等于 $n$ 的整数 $k$ , $\omega _ {n} ^ {k}=\omega _ {n} ^ {k\ mod\ n}$ 。

单位复数根的性质

  接下来介绍单位复数根的几个性质:
  消去引理:对任何整数 $n≥0,\ k≥0,\ d>0$ ,存在:

$$\omega _ {dn} ^ {dk}=\omega _ {n} ^ {k}$$

  这个道理是比较明显的。

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

  折半引理:如果 $n>0$ 且为偶数,那么 $n$ 个 $n$ 次单位复数根的平方的集合就是 $\frac{n}{2}$ 个 $\frac{n}{2}$ 次单位复数根的集合。
  根据消去引理,对任意非负整数 $k$ ,有 $(\omega _ {n} ^ {k})^2=\omega _ {n/2} ^ {k}$ 。如果对所有 $n$ 次单位复数根进行平方,那么获得每个 $n/2$ 次单位根正好 $2$ 次,因为 $\omega _ {n} ^ {k}$ 与 $\omega _ {n} ^ {k+n/2}$ 的平方相同。

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

  这个理论对接下来的计算非常重要。

离散傅里叶变换DFT

  回到正题。我们希望计算次数界为 $n$ 的多项式 $A(x)=\sum\limits_{j=0}^{n-1}{a_jx^j}$ 在 $n$ 个 $n$ 次单位复数根处的取值。 $A(x)$ 以系数形式给出 $a=(a_0,\ a_1,\ \cdots,\ a_{n-1})$ 。对每一个单位根对应的结果,我们知道结果 $y_k$ :

$$y_k=A(\omega _ {n} ^ {k})=\sum\limits_{j=0}^{n-1}{a_j\omega _ n ^{kj}}$$

  得到的向量 $y=(y_0,\ y_1,\ \cdots,\ y_{n-1})$ 称为系数向量 $a=(a_0,\ a_1,\ \cdots,\ a_{n-1})$ 的离散傅里叶变换(DFT)

快速傅里叶变换FFT

  从下文开始我们假设 $n$ 是 $2$ 的整数幂。
  对于一个多项式 $A(x)=\sum\limits_{j=0}^{n-1}{a_jx^j}$ ,我们可以按照系数的奇偶性对他进行处理:
\begin{align}
A(x)&=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1}\\
&=(a_0+a_2x^2+a_4x^4+\cdots +a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdots +a_{n-1}x^{n-1})\\
&=(a_0+a_2x^2+a_4x^4+\cdots +a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots +a_{n-1}x^{n-2})\\
\end{align}

  发现有些地方有些相似?我们再处理一下。假设两个多项式:
\begin{align}
A^{[0]}(x)&=(a_0+a_2x+a_4x^2+\cdots +a_{n-2}x^{n/2-1})\\
A^{[1]}(x)&=(a_1+a_3x+a_5x^2+\cdots +a_{n-1}x^{n/2-1})
\end{align}

  显然有如下结论:

$$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$$

  现在问题转化为:求出次数界为 $n/2$ 的多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$在 $(\omega _ {n} ^ {0})^2,\ (\omega _ {n} ^ {1})^2,\ \cdots,\ (\omega _ {n} ^ {n-1})^2$ 的取值。
  然后我们利用折半引理,发现上述 $n$ 个式子的值其实仅由 $n/2$ 个 $n/2$ 次单位根组成,每个根出现两次。
  然后之后的求值便是一个递归问题了。代码大概如下:

cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
void fft(complex<double> *a,int n){
//c++支持附属模板complex
if (n==1) return; //如果就一个常数项就返回
int mid=n>>1;
complex<double> tmp[maxn];
//tmp作为临时数组方便a的处理
for (int i=0;i<mid;++i){
tmp[i]=a[i*2]; tmp[mid+i]=a[i*2+1];
//偶数项在0~mid-1,奇数在mid~n-1
}
for (int i=0;i<n;++i) a[i]=tmp[i];
//分治
fft(a,mid); fft(a+mid,mid);
//以下部分将在下面讲述
for (int i=0;i<mid;++i){
complex<double> omega=(cos(2*pi*i/n),sin(2*pi*i/n));
tmp[i]=a[i]+omega*a[i]; tmp[i+mid]=a[i]-omega*a[i+mid];
}
for (int i=0;i<n;++i) a[i]=tmp[i];
}

  关于第 $17$ 行的两个操作:

$$y_k=A(\omega _ {n} ^ {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_{k+(n/2)}=A(\omega _ {n} ^ {k+(n/2)}) = A^{[0]}(\omega _ {n} ^ {2k+n})+\omega _ n ^ {k+(n/2)} A^{[1]} (\omega _ {n} ^ {2k+n})=A^{[0]}(\omega _ {n/2} ^ {k})+\omega _ n ^ {k+(n/2)} A^{[1]} (\omega _ {n/2} ^ {k})$$

由于:

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

得到:

$$y_{k+(n/2)}=A(\omega _ {n} ^ {k+(n/2)}) =A^{[0]}(\omega _ {n/2} ^ {k})-\omega _ n ^ k A^{[1]} (\omega _ {n/2} ^ {k})$$

快速傅里叶逆变换IFFT

  处理完之后,我们需要把点值表达的多项式再转为系数表达的。所以我们需要再在单位复数根处插值,即由点值表达转变为系数表达。
  我们可以把 DFT 写成一个矩阵方程,观察其逆矩阵的形式: $y=V_n a$ 。
\begin{equation}
\left[\begin{array}{a1}
y_0\\
y_1\\
y_2\\
y_3\\
\vdots\\
y_{n-1}\\
\end{array}\right]
=
\left[\begin{array}{a1}
1 &1 &1 &1 &\cdots &1\\
1 &\omega _ {n} ^ {1} &\omega _ {n} ^ {2} &\omega _ {n} ^ {3} &\cdots &\omega _ {n} ^ {n-1}\\
1 &\omega _ {n} ^ {2} &\omega _ {n} ^ {4} &\omega _ {n} ^ {6} &\cdots &\omega _ {n} ^ {2(n-1)}\\
1 &\omega _ {n} ^ {3} &\omega _ {n} ^ {6} &\omega _ {n} ^ {9} &\cdots &\omega _ {n} ^ {3(n-1)}\\
\vdots &\vdots &\vdots &\vdots &\ddots &\vdots\\
1 &\omega _ {n} ^ {n-1} &\omega _ {n} ^ {2(n-1)} &\omega _ {n} ^ {3(n-1)} &\cdots &\omega _ {n} ^ {(n-1)(n-1)}\\
\end{array}\right]
\left[\begin{array}{a1}
a_0\\
a_1\\
a_2\\
a_3\\
\vdots\\
a_{n-1}\\
\end{array}\right]
\end{equation}

  对于逆运算,我们通过 $y$ 乘以 $V_n$ 的逆矩阵 $V_n^{-1}$ 进行处理。由表格我们可以清楚地看出:对于 $j,\ k=0,\ 1,\ 2,\ \cdots,\ n-1$ , $V_n$ 的 $(k,\ j)$ 处元素为 $\omega _ n ^ {kj}$ 。
  接下来送命题:矩阵求逆。
  通过奇妙的手段,我们可以得到: $V_n^{-1}$ 的 $(j,\ k)$ 处的元素是 $\omega _ n ^ {-kj}/n$ 。得到最后系数表达的向量 $a_j=\frac{1}{n}\sum\limits_{k=0}^{n-1}{y_k\omega _ n ^{-kj}}$ 。
  有没有发现 $a_j=\frac{1}{n}\sum\limits_{k=0}^{n-1}{y_k\omega _ n ^{-kj}}$ 和 $y_k=\sum\limits_{j=0}^{n-1}{a_j\omega _ n ^{kj}}$ 比较相似?抛开前面不管,最后的单位根其实互为倒数(感觉是废话)。
  我们注意到单位根有一个非常好的性质:实部为 $a$ ,虚部为 $b$ ,则 $a^2+b^2=1$ ($cos^2x+sin^2x=1$)。
  我们还知道,一个复数 $z$ 与他的共轭复数 $\bar z$ 相乘,得到的结果就是 $a^2+b^2=1$ ,与我们想要找的倒数契合。
  所以 $a_j=\sum\limits_{k=0}^{n-1}{y_k\omega _ n ^{-kj}}$ 是很方便求出的:原本乘单位复数根的地方乘以原来单位根的共轭复数即可
  还有一个 $\frac{1}{n}$ ?最后乘上就行。

代码(simple)

  由上文可以得出 FFT 与 IFFT 是可以写在一块的。我们可以这样写:

FFT函数:

cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void fft(complex<double> *a,int n,int inv){
//inv代表是FFT还是IFFT(即是否要共轭),值为1或-1
if (n==1) return;
int mid=n>>1;
complex<double> tmp[maxn];
for (int i=0;i<mid;++i){
tmp[i]=a[i*2]; tmp[mid+i]=a[i*2+1];
}
for (int i=0;i<n;++i) a[i]=tmp[i];
fft(a,mid,inv); fft(a+mid,mid,inv);
for (int i=0;i<mid;++i){
complex<double> omega=(cos(2*pi*i/n),inv*sin(2*pi*i/n));
//乘上-1即代表求共轭
tmp[i]=a[i]+omega*a[i]; tmp[i+mid]=a[i]-omega*a[i+mid];
}
for (int i=0;i<n;++i) a[i]=tmp[i];
}

主要步骤:

cpp
1
2
3
4
5
6
7
8
9
complex<double> a[maxn],b[maxn];
int ans[maxn];
void work(){
fft(a,n,1); fft(b,n,1); //FFT转点值
for (int i=0;i<n;++i) a[i]*=b[i]; //点值相乘
fft(a,n,-1); //IFFT转系数
for (int i=0;i<n;++i) c[i]=(a[i].real()/n+0.5);
//别忘了最后的1/n,同时注意精度
}

FFT的优化

  其实我们目前得到的FFT是不算特别健全的。下面讨论两种优化方法。

迭代FFT

  一般来说递归比迭代的效率是要低的(而且中间开个大数组)。我们就考虑不用递归。
  我们画一个表格来模拟一下系数向量:

$ a_0$,$\ a_1$,$\ a_2$,$\ a_3$,$\ a_4$,$\ a_5$,$\ a_6$,$\ a_7$
$a_0$,$\ a_2$,$\ a_4$,$\ a_6$ $a_1$,$\ a_3$,$\ a_5$,$\ a_7$
$a_0$,$\ a_4$ $a_2$,$\ a_6$ $a_1$,$\ a_5$ $a_3$,$\ a_7$
$a_0$ $a_4$ $a_2$ $a_6$ $a_1$ $a_5$ $a_3$ $a_7$

  没有发现?我们再来转变一下(为了偷懒蒟蒻博主省去a了只写下标的二进制):

$000$,$\ 001$,$\ 010$,$\ 011$,$\ 100$,$\ 101$,$\ 110$,$\ 111$
$000$,$\ 010$,$\ 100$,$\ 110$ $001$,$\ 011$,$\ 101$,$\ 111$
$000$,$\ 100$ $010$,$\ 110$ $001$,$\ 101$ $011$,$\ 111$
$000$ $100$ $010$ $110$ $001$ $101$ $011$ $111$

  我们可以发现 $i$ 号元素最后的位置在 $i$ 的二进制倒序的结果。
  所以我们可以先找到对应的位置,然后再往上进行还原就ok了。

蝴蝶操作

  之前代码前我们开了个数组,主要是为了下面代码

cpp
1
2
3
4
5
for (int i=0;i<mid;++i){
complex<double> omega=(cos(2*pi*i/n),inv*sin(2*pi*i/n));
tmp[i]=a[i]+omega*a[i]; tmp[i+mid]=a[i]-omega*a[i+mid];
}
for (int i=0;i<n;++i) a[i]=tmp[i];

  使用这个临时数组防止 $a[i]$ 的值被覆盖了。那我们换个顺序,或者提前存一下不就行了?

cpp
1
2
3
4
5
for (int i=0;i<mid;++i){
complex<double> omega=(cos(2*pi*i/n),inv*sin(2*pi*i/n));
a[i+mid]=a[i]-omega*a[i+mid];
a[i]=a[i]+omega*a[i];
}

蒟蒻代码

  以洛谷的这道模板题开刀吧。P3803 【模板】多项式乘法(FFT)

cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
#define ll long long
#define il inline
#define db double
#define cp complex_number
const double PI=acos(-1.0);
const int maxn=5000000;
il int read(){
int x=0,f=1; char ch=getchar();
while (ch<'0'||ch>'9'){
if (ch=='-') f=-1; ch=getchar();
}
while (ch>='0'&&ch<='9'){
x=(x<<1)+(x<<3)+ch-'0'; ch=getchar();
}
return x*f;
}
struct complex_number{
db rel,img;
complex_number(double x=0,double y=0){
rel=x; img=y;
}
};
cp operator + (cp a,cp b){
return cp(a.rel+b.rel,a.img+b.img);
}
cp operator - (cp a,cp b){
return cp(a.rel-b.rel,a.img-b.img);
}
cp operator * (cp a,cp b){
return cp(a.rel*b.rel-a.img*b.img,a.rel*b.img+b.rel*a.img);
}
cp aa[maxn],bb[maxn];
int n,m,limits=1,lim=0;
int rev[maxn];
void fft(cp *a,int len,int inv){
for (int i=0;i<len;++i){
if (i<rev[i]) swap(a[i],a[rev[i]]);
}
for (int mid=1;mid<len;mid<<=1){
cp tmp=cp(cos(PI/mid),inv*sin(PI/mid));
for (int i=0;i<len;i+=mid<<1){
cp omega=cp(1,0);
for (int j=0;j<mid;++j,omega=omega*tmp){
cp x=a[i+j],y=omega*a[i+j+mid];
a[i+j]=x+y; a[i+j+mid]=x-y;
}
}
}
}
int main(){
n=read(); m=read();
while (limits<=m+n) limits<<=1,++lim;
for (int i=0;i<=n;++i) aa[i].rel=read();
for (int i=0;i<=m;++i) bb[i].rel=read();
for (int i=0;i<limits;++i){
rev[i]=(rev[i>>1]>>1)|((i&1)<<(lim-1));
}
fft(aa,limits,1);
fft(bb,limits,1);
for (int i=0;i<=limits;++i){
aa[i]=aa[i]*bb[i];
}
fft(aa,limits,-1);
for (int i=0;i<=n+m;++i){
printf("%d ",(int)(aa[i].rel/limits+0.5));
}
return 0;
}

参考资料

文章作者: lonlyn
文章链接: https://lonlyn.gitee.io/blog/2019/11/21/study_fft/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 lonlyn's blog

评论