概述
多项式乘法是一种很常见的运算。朴素算法的时间复杂度为
O
(
n
2
)
O(n^2)
O(n2),Karatsuba 乘法的时间复杂度为
O
(
n
log
2
3
)
O(n^{{\log_2}3})
O(nlog23)但是常数爆炸。下面介绍的快速傅里叶变换(FFT)可以在
O
(
n
log
n
)
O(n\log n)
O(nlogn)的时间内计算两个多项式的乘法。
注:以下多项式均为一元多项式
前置技能
复数
记
i
2
=
−
1
i^2=-1
i2=−1,
i
i
i为虚数单位。
形如
a
+
b
i
(
a
,
b
∈
R
)
a+bi(a,b \in \mathbb{R})
a+bi(a,b∈R)的数称为复数,其中
a
a
a称为实部,
b
b
b称为虚部。
a
2
+
b
2
\sqrt{a^2+b^2}
a2+b2称为
z
z
z的模长,记作
∣
z
∣
|z|
∣z∣。
复平面
建立平面直角坐标系
x
O
y
xOy
xOy,对于复数
z
=
a
+
b
i
z=a+bi
z=a+bi,点
Z
(
a
,
b
)
Z(a,b)
Z(a,b)为其在坐标系内对应的点,
O
Z
→
=
(
a
,
b
)
\overrightarrow{OZ}=(a,b)
OZ=(a,b)为其对应的向量。
显然一定存在
θ
\theta
θ,使
z
=
cos
θ
+
i
sin
θ
z=\cos\theta+i\sin\theta
z=cosθ+isinθ。
θ
\theta
θ称为
z
z
z的辐角,记作
A
r
g
(
z
)
Arg(z)
Arg(z),当
θ
∈
(
−
π
,
π
)
\theta\in(-\pi,\pi)
θ∈(−π,π)时,
θ
\theta
θ称为
z
z
z的辐角主值,记作
a
r
g
(
z
)
arg(z)
arg(z)。
复数加法
对于
z
1
=
a
1
+
b
1
i
,
z
2
=
a
2
+
b
2
i
z_1=a_1+b_1i,z_2=a_2+b_2i
z1=a1+b1i,z2=a2+b2i,有
z
1
+
z
2
=
(
a
1
+
a
2
)
+
(
b
1
+
b
2
)
i
z_1+z_2=(a_1+a_2)+(b_1+b_2)i
z1+z2=(a1+a2)+(b1+b2)i。
即在复平面内对应的向量相加。
复数减法
对于
z
1
=
a
1
+
b
1
i
,
z
2
=
a
2
+
b
2
i
z_1=a_1+b_1i,z_2=a_2+b_2i
z1=a1+b1i,z2=a2+b2i,有
z
1
−
z
2
=
(
a
1
−
a
2
)
+
(
b
1
−
b
2
)
i
z_1-z_2=(a_1-a_2)+(b_1-b_2)i
z1−z2=(a1−a2)+(b1−b2)i。
即在复平面内对应的向量相减。
复数乘法
对于
z
1
=
a
1
+
b
1
i
,
z
2
=
a
2
+
b
2
i
z_1=a_1+b_1i,z_2=a_2+b_2i
z1=a1+b1i,z2=a2+b2i,有
z
1
z
2
=
(
a
1
a
2
−
b
1
b
2
)
+
(
a
1
b
2
+
a
2
b
2
)
i
z_1z_2=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_2)i
z1z2=(a1a2−b1b2)+(a1b2+a2b2)i。
可以证明,复数乘法的几何意义为模长相乘,辐角相加,如下图所示。
单位根
方程
z
n
=
1
z^n=1
zn=1的复数根称为
n
n
n次单位根。
n
n
n次单位根有
n
n
n个,
cos
2
k
π
n
+
i
sin
2
k
π
n
(
0
≤
k
≤
n
−
1
,
k
∈
Z
)
\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}(0\leq k\leq n-1,k\in\mathbb{Z})
cosn2kπ+isinn2kπ(0≤k≤n−1,k∈Z)。
记
ω
n
=
cos
2
π
n
+
i
sin
2
π
n
\omega_n=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}
ωn=cosn2π+isinn2π,则
n
n
n次单位根就是
ω
n
0
,
ω
n
1
,
ω
n
2
,
⋯
 
,
ω
n
n
−
1
\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}
ωn0,ωn1,ωn2,⋯,ωnn−1。
证明:
显然有
∣
ω
n
k
∣
=
cos
2
2
k
π
n
+
sin
2
2
k
π
n
=
1
|\omega_n^k|=\sqrt{\cos^2\frac{2k\pi}{n}+\sin^2\frac{2k\pi}{n}}=1
∣ωnk∣=cos2n2kπ+sin2n2kπ=1。考虑复数乘法的几何意义,
(
ω
n
k
)
n
{(\omega_n^k)}^n
(ωnk)n的模长为
1
1
1,辐角为
2
k
π
n
⋅
n
=
2
k
π
\frac{2k\pi}{n}\cdot n=2k\pi
n2kπ⋅n=2kπ,故
(
ω
n
k
)
n
=
cos
2
k
π
+
i
sin
2
k
π
=
1
{(\omega_n^k)}^n=\cos2k\pi+i\sin2k\pi=1
(ωnk)n=cos2kπ+isin2kπ=1。
单位根的性质
- ω n k = ω d n d k \omega_n^k=\omega_{dn}^{dk} ωnk=ωdndk
- ω n k = − ω n k + n 2 \omega_n^k=-\omega_n^{k+\frac{n}{2}} ωnk=−ωnk+2n
- ω n 0 = ω n n = 1 \omega_n^0=\omega_n^n=1 ωn0=ωnn=1
证明:
还是考虑复数乘法的几何意义。
- ω n k \omega_n^k ωnk与 ω d n d k \omega_{dn}^{dk} ωdndk在复平面内对应的向量相等,且模长都为 1 1 1,故 ω n k = ω d n d k \omega_n^k=\omega_{dn}^{dk} ωnk=ωdndk。
- ω n k \omega_n^k ωnk与 ω n k + n 2 \omega_n^{k+\frac{n}{2}} ωnk+2n在复平面内对应的向量关于坐标原点对称,故 ω n k = − ω n k + n 2 \omega_n^k=-\omega_n^{k+\frac{n}{2}} ωnk=−ωnk+2n。
- 显然成立。
复数的代码实现
struct Complex {
double a, b;
Complex() { a = b = 0; }
Complex(double a, double b) : a(a), b(b) {}
};
Complex operator+(Complex x, Complex y) { return Complex(x.a + y.a, x.b + y.b); }
Complex operator-(Complex x, Complex y) { return Complex(x.a - y.a, x.b - y.b); }
Complex operator*(Complex x, Complex y) { return Complex(x.a*y.a - x.b*y.b, x.a*y.b + x.b*y.a); }
多项式的系数表示法与点值表示法
系数表示法
通常,形如 A ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n A(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n A(x)=a0+a1x+a2x2+⋯+anxn叫做多项式的点值表示法,即最常见的表示法。该多项式的次数为 n n n,通常记作 d e g A = n degA=n degA=n。
点值表示法
如果选取
n
+
1
n+1
n+1个数带入多项式进行求值,得到
A
(
x
0
)
,
A
(
x
1
)
,
⋯
 
,
A
(
x
n
)
A(x_0),A(x_1),\cdots,A(x_n)
A(x0),A(x1),⋯,A(xn),则称
A
(
x
)
=
{
(
x
i
,
A
(
x
i
)
)
∣
0
≤
i
≤
n
,
i
∈
Z
}
A(x)=\{(x_i,A(x_i))|0 \leq i \leq n ,i \in \mathbb{Z}\}
A(x)={(xi,A(xi))∣0≤i≤n,i∈Z}为多项式
A
(
x
)
A(x)
A(x)的点值表示法。
点值表示法的优越性在于可以在
O
(
n
)
O(n)
O(n)的时间内计算多项式乘法。
例如:计算
A
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
⋯
+
a
n
x
n
A(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n
A(x)=a0+a1x+a2x2+⋯+anxn与
B
(
x
)
=
b
0
+
b
1
x
+
b
2
x
2
+
⋯
+
b
n
x
n
B(x)=b_0+b_1x+b_2x^2+\cdots+b_nx^n
B(x)=b0+b1x+b2x2+⋯+bnxn的乘积
C
(
x
)
C(x)
C(x)。
在系数表示法中,
C
(
x
)
=
A
(
x
)
B
(
x
)
=
∑
i
=
0
2
n
c
i
x
i
C(x)=A(x)B(x)=\sum\limits_{i=0}^{2n}c_ix^i
C(x)=A(x)B(x)=i=0∑2ncixi,其中
c
i
=
∑
j
=
0
i
a
j
b
i
−
j
c_i=\sum\limits_{j=0}^{i}a_jb_{i-j}
ci=j=0∑iajbi−j。不难看出,总计算次数为
O
(
n
2
)
O(n^2)
O(n2)。
而在系数表示法中,
C
(
x
)
=
{
(
x
i
,
A
(
x
i
)
B
(
x
i
)
)
∣
0
≤
i
≤
n
,
i
∈
Z
}
C(x)=\{(x_i,A(x_i)B(x_i))|0 \leq i \leq n ,i \in \mathbb{Z}\}
C(x)={(xi,A(xi)B(xi))∣0≤i≤n,i∈Z},总计算次数为
O
(
n
)
O(n)
O(n)。
快速傅里叶变换
前面已经提到过,点值表示法可以在 O ( n ) O(n) O(n)的时间内计算两个多项式的乘积,所以时间复杂度的瓶颈在于点值表示法与系数表示法直接的转化,即多项式的求值与插值。快速傅里叶变换可以在 O ( n log n ) O(n\log n) O(nlogn)的时间完成求值与插值操作。
求值
现在的问题是,如何选取代入的
n
n
n个点,才能实现快速求值。
考虑到单位根的性质,将
n
n
n次单位根的
n
n
n个值代入多项式,可以简化计算。
设
n
−
1
n-1
n−1次多项式
A
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
⋯
a
n
−
1
x
n
−
1
A(x)=a_0+a_1x+a_2x^2+\cdots a_{n-1}x^{n-1}
A(x)=a0+a1x+a2x2+⋯an−1xn−1。
n
=
2
m
(
m
∈
Z
)
n=2^m(m\in\mathbb{Z})
n=2m(m∈Z),如果不满足则将
n
n
n补成
2
m
2^m
2m的形式。
将
A
(
x
)
A(x)
A(x)按照
x
x
x次数的奇偶性分类。
设
A
0
(
x
)
=
a
0
+
a
2
x
+
a
4
x
2
+
⋯
+
a
n
−
2
x
n
2
−
1
A_0(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1}
A0(x)=a0+a2x+a4x2+⋯+an−2x2n−1
,
A
1
(
x
)
=
a
1
+
a
3
x
+
a
5
x
2
+
⋯
+
a
n
−
1
x
n
2
−
1
A_1(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1}
A1(x)=a1+a3x+a5x2+⋯+an−1x2n−1。
可以看出
A
(
x
)
=
A
0
(
x
2
)
+
x
A
1
(
x
2
)
A(x)=A_0(x^2)+xA_1(x^2)
A(x)=A0(x2)+xA1(x2)。
现在,将
ω
n
k
(
k
<
n
2
)
\omega_n^k(k<\frac{n}{2})
ωnk(k<2n)代入,得
A
(
ω
n
k
)
=
A
0
(
ω
n
2
k
)
+
ω
n
k
A
1
(
ω
n
2
k
)
A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^kA_1(\omega_n^{2k})
A(ωnk)=A0(ωn2k)+ωnkA1(ωn2k)
=
A
0
(
ω
n
2
k
)
+
ω
n
k
A
1
(
ω
n
2
k
)
= A_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k})
=A0(ω2nk)+ωnkA1(ω2nk)
同理,将
ω
n
k
+
n
2
\omega_n^{k+\frac{n}{2}}
ωnk+2n代入,得
A
(
ω
n
k
+
n
2
)
=
A
0
(
ω
n
2
k
+
n
)
+
ω
n
k
+
n
2
A
1
(
ω
n
2
k
+
n
)
A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_1(\omega_n^{2k+n})
A(ωnk+2n)=A0(ωn2k+n)+ωnk+2nA1(ωn2k+n)
=
A
0
(
ω
n
2
k
+
n
2
)
−
ω
n
k
A
1
(
ω
n
2
k
+
n
2
)
=A_0(\omega_{\frac{n}{2}}^{k+\frac{n}{2}})-\omega_n^kA_1(\omega_{\frac{n}{2}}^{k+\frac{n}{2}})
=A0(ω2nk+2n)−ωnkA1(ω2nk+2n)
=
A
0
(
ω
n
2
k
)
−
ω
n
k
A
1
(
ω
n
2
k
)
=A_0(\omega_{\frac{n}{2}}^k)-\omega_n^kA_1(\omega_{\frac{n}{2}}^k)
=A0(ω2nk)−ωnkA1(ω2nk)
可以发现两个式子化简后的结果只有符号的区别,且可以划分为两个子问题解决,于是可以考虑分治(代码暂时不贴)。时间复杂度为 O ( n log n ) O(n\log n) O(nlogn)
插值
现在要把点值表达式还原成系数表达式。若使用拉格朗日插值,时间复杂度为
O
(
n
2
)
O(n^2)
O(n2),需要优化。
考虑这是一个解方程的问题。
{
a
0
(
ω
n
0
)
0
+
a
1
(
ω
n
0
)
1
+
⋯
+
a
n
−
1
(
ω
n
0
)
n
−
1
=
A
(
ω
n
0
)
a
0
(
ω
n
1
)
0
+
a
1
(
ω
n
1
)
1
+
⋯
+
a
n
−
1
(
ω
n
1
)
n
−
1
=
A
(
ω
n
1
)
⋮
⋮
⋮
⋮
a
0
(
ω
n
n
−
1
)
0
+
a
1
(
ω
n
n
−
1
)
1
+
⋯
+
a
n
−
1
(
ω
n
n
−
1
)
n
−
1
=
A
(
ω
n
n
−
1
)
\begin{cases} 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)\\ \vdots & \vdots &\vdots&\vdots\\ 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{cases}
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧a0(ωn0)0+a0(ωn1)0+⋮a0(ωnn−1)0+a1(ωn0)1+a1(ωn1)1+⋮a1(ωnn−1)1+⋯⋯⋮⋯+an−1(ωn0)n−1+an−1(ωn1)n−1⋮+an−1(ωnn−1)n−1=A(ωn0)=A(ωn1)=A(ωnn−1)
写成矩阵的形式
[
(
ω
n
0
)
0
(
ω
n
0
)
1
⋯
(
ω
n
0
)
n
−
1
(
ω
n
1
)
0
(
ω
n
1
)
1
⋯
(
ω
n
1
)
n
−
1
⋮
⋮
⋱
⋮
(
ω
n
n
−
1
)
0
(
ω
n
n
−
1
)
1
⋯
(
ω
n
n
−
1
)
n
−
1
]
[
a
0
a
1
⋮
a
n
−
1
]
=
[
A
(
ω
n
0
)
A
(
ω
n
1
)
⋮
A
(
ω
n
n
−
1
)
]
\begin{bmatrix} {(\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}\\ \vdots & \vdots & \ddots & \vdots\\ {(\omega_n^{n-1})}^0 & {(\omega_n^{n-1})}^1& \cdots & {(\omega_n^{n-1})}^{n-1}\\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \\ \end{bmatrix}= \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \\ \end{bmatrix}
⎣⎢⎢⎢⎢⎡(ωn0)0(ωn1)0⋮(ωnn−1)0(ωn0)1(ωn1)1⋮(ωnn−1)1⋯⋯⋱⋯(ωn0)n−1(ωn1)n−1⋮(ωnn−1)n−1⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡A(ωn0)A(ωn1)⋮A(ωnn−1)⎦⎥⎥⎥⎤
现在可以使用高斯消元,时间复杂度为
O
(
n
3
)
O(n^3)
O(n3),完美的负优化。。。
尝试去直接构造结果,设系数矩阵为
V
V
V,考虑这样一个矩阵
D
=
[
(
ω
n
−
0
)
0
(
ω
n
−
0
)
1
⋯
(
ω
n
−
0
)
n
−
1
(
ω
n
−
1
)
0
(
ω
n
−
1
)
1
⋯
(
ω
n
−
1
)
n
−
1
⋮
⋮
⋱
⋮
(
ω
n
−
(
n
−
1
)
)
0
(
ω
n
−
(
n
−
1
)
)
1
⋯
(
ω
n
−
(
n
−
1
)
)
n
−
1
]
D= \begin{bmatrix} {(\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}\\ \vdots & \vdots & \ddots & \vdots\\ {(\omega_n^{-(n-1)})}^0 & {(\omega_n^{-(n-1)})}^1& \cdots & {(\omega_n^{-(n-1)})}^{n-1}\\ \end{bmatrix}
D=⎣⎢⎢⎢⎢⎡(ωn−0)0(ωn−1)0⋮(ωn−(n−1))0(ωn−0)1(ωn−1)1⋮(ωn−(n−1))1⋯⋯⋱⋯(ωn−0)n−1(ωn−1)n−1⋮(ωn−(n−1))n−1⎦⎥⎥⎥⎥⎤
设
E
=
D
V
E=DV
E=DV,则
e
i
j
=
∑
k
=
0
n
−
1
d
i
k
v
k
j
=
∑
k
=
0
n
−
1
ω
n
−
i
k
ω
n
k
j
=
∑
k
=
0
n
−
1
ω
n
k
(
j
−
i
)
e_{ij}=\sum\limits_{k=0}^{n-1}d_{ik}v_{kj}=\sum\limits_{k=0}^{n-1}\omega_n^{-ik}\omega_n^{kj}=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}
eij=k=0∑n−1dikvkj=k=0∑n−1ωn−ikωnkj=k=0∑n−1ωnk(j−i)
当
j
=
i
j=i
j=i时,
e
i
j
=
n
e_{ij}=n
eij=n;
当
j
̸
=
i
j\not=i
j̸=i时,
e
i
j
=
∑
k
=
0
n
−
1
ω
n
k
(
j
−
i
)
=
1
−
(
ω
n
j
−
i
)
n
1
−
ω
n
j
−
i
=
0
e_{ij}=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}=\frac{1-{(\omega_n^{j-i})}^n}{1-\omega_n^{j-i}}=0
eij=k=0∑n−1ωnk(j−i)=1−ωnj−i1−(ωnj−i)n=0。
故
I
n
=
1
n
E
I_n=\frac{1}{n}E
In=n1E,所以
1
n
D
=
V
−
1
\frac{1}{n}D=V^{-1}
n1D=V−1。
此时可以得到
[
a
0
a
1
⋮
a
n
−
1
]
=
1
n
[
(
ω
n
−
0
)
0
(
ω
n
−
0
)
1
⋯
(
ω
n
−
0
)
n
−
1
(
ω
n
−
1
)
0
(
ω
n
−
1
)
1
⋯
(
ω
n
−
1
)
n
−
1
⋮
⋮
⋱
⋮
(
ω
n
−
(
n
−
1
)
)
0
(
ω
n
−
(
n
−
1
)
)
1
⋯
(
ω
n
−
(
n
−
1
)
)
n
−
1
]
[
A
(
ω
n
0
)
A
(
ω
n
1
)
⋮
A
(
ω
n
n
−
1
)
]
\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \\ \end{bmatrix}=\frac{1}{n} \begin{bmatrix} {(\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}\\ \vdots & \vdots & \ddots & \vdots\\ {(\omega_n^{-(n-1)})}^0 & {(\omega_n^{-(n-1)})}^1& \cdots & {(\omega_n^{-(n-1)})}^{n-1}\\ \end{bmatrix} \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \\ \end{bmatrix}
⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤=n1⎣⎢⎢⎢⎢⎡(ωn−0)0(ωn−1)0⋮(ωn−(n−1))0(ωn−0)1(ωn−1)1⋮(ωn−(n−1))1⋯⋯⋱⋯(ωn−0)n−1(ωn−1)n−1⋮(ωn−(n−1))n−1⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡A(ωn0)A(ωn1)⋮A(ωnn−1)⎦⎥⎥⎥⎤
所以,插值等价于将单位根的倒数代入求值的结果除以
n
n
n。
这样,求值和插值都可以在
O
(
n
log
n
)
O(n\log n)
O(nlogn)的时间内完成了。