这是一篇基于OI Wiki中的内容进行整理的算法竞赛中关于多项式的相关内容的第一篇笔记。
1.基本概念
对于求和式 ∑ a n x n \displaystyle\sum\limits a_nx^n ∑anxn，如果是有限项相加，则这个式子被称作多项式(Polynomial)，记作 f ( x ) = ∑ i = 0 n a n x n f(x)=\displaystyle\sum\limits_{i=0}^n a_nx^n f(x)=i=0∑nanxn。
多项式还可以写成这样的形式： f ( x ) = f 0 + f 1 x + f 2 x 2 + ⋯ + f n x n f(x)=f_0+f_1x+f_2x^2+\cdots +f_nx^n f(x)=f0+f1x+f2x2+⋯+fnxn。
在这个 f ( x ) f(x) f(x)中，可以认为每一项的系数 f 0 , f 1 , ⋯ , f n f_0,f_1,\cdots, f_n f0,f1,⋯,fn是这个多项式中的变量，不同的系数表示了不同的多项式，即 f 0 , f 1 , ⋯ , f n f_0,f_1,\cdots, f_n f0,f1,⋯,fn是 f f f的特征。如果将一个多项式当作一个元素，那么也可以写成这样的形式： f = ⟨ f 0 , f 1 , f 2 , ⋯ f n ⟩ f=\langle f_0,f_1,f_2,\cdots f_n\rangle f=⟨f0,f1,f2,⋯fn⟩，如果多项式 f f f是 R \mathbb{R} R上的多项式，那么 f 0 , f 1 , ⋯ , f n ∈ R f_0,f_1,\cdots,f_n∈\mathbb{R} f0,f1,⋯,fn∈R。
对于一个多项式
f
(
x
)
f(x)
f(x)，称其最高次项的次数为该多项式的度(Degree)，也称作次数。记作
deg
f
\deg f
degf。
例如，
f
(
x
)
=
3
x
2
+
6
x
+
7
f(x)=3x^2+6x+7
f(x)=3x2+6x+7，这个多项式的最高次项是
2
2
2，则称这是一个
2
2
2度的多项式。
多项式乘法
关于多项式，最核心的操作就是两个多项式的乘法，即将 R \mathbb{R} R上的两个多项式 f ( x ) f(x) f(x)和 g ( x ) g(x) g(x)相乘。
若有
f
(
x
)
=
a
0
+
a
1
x
+
⋯
+
a
n
x
n
f(x)=a_0+a_1x+\cdots+a_nx^n
f(x)=a0+a1x+⋯+anxn，
g
(
x
)
=
b
0
+
b
1
x
+
⋯
+
b
n
x
m
g(x)=b_0+b_1x+\cdots+b_nx^m
g(x)=b0+b1x+⋯+bnxm，那么我们的目标就是求出
R
\mathbb{R}
R上的多项式
Q
(
x
)
=
f
(
x
)
⋅
g
(
x
)
Q(x)=f(x)\cdot g(x)
Q(x)=f(x)⋅g(x)：
Q
(
x
)
=
∑
i
=
0
n
∑
j
=
0
m
a
i
b
j
x
i
+
j
=
c
0
+
c
1
x
+
⋯
+
c
n
+
m
x
n
+
m
Q(x)=\displaystyle\sum\limits_{i=0}^n \displaystyle\sum\limits_{j=0}^m a_ib_jx^{i+j}=c_0+c_1x+\cdots+c_{n+m}x^{n+m}
Q(x)=i=0∑nj=0∑maibjxi+j=c0+c1x+⋯+cn+mxn+m
显然，逐项相乘求得 Q ( x ) Q(x) Q(x)各项系数的方法的时间复杂度是 O ( n m ) O(nm) O(nm)。一般情况下，这个时间复杂度我们是不能接受的。接下来的内容，就是利用各种手段，将求解多项式操作的时间复杂度变得可以接受。
复数
定义 i 2 = − 1 i^2=-1 i2=−1。若 z z z是复数域 C \mathbb{C} C上的复数，那么可以写作： z = x + y i , x , y ∈ R z=x+yi ,x,y∈\mathbb{R} z=x+yi,x,y∈R。
我们可以用一个 x x x轴表示 R \mathbb{R} R上的所有实数，同理我们也可以用平面直角坐标系表示 C \mathbb{C} C上的一个复数，对应关系是：平面直角坐标系（我们称为复平面）上 ( x , y ) (x,y) (x,y)点表示复数 z = x + y i z=x+yi z=x+yi。
如图所示，连接
(
0
,
0
)
(0,0)
(0,0)和
(
x
,
y
)
(x,y)
(x,y)，构成一个向量，设这个向量与
x
x
x轴的夹角为
θ
\theta
θ。
复数的模
∣
z
∣
=
x
2
+
y
2
\lvert z\rvert=\sqrt{x^2+y^2}
∣z∣=x2+y2，坐标系上就是
(
0
,
0
)
(0,0)
(0,0)到
(
x
,
y
)
(x,y)
(x,y)的距离。
复数还有以下形式：
z
=
x
+
y
i
=
∣
z
∣
(
cos
θ
+
i
sin
θ
)
=
∣
z
∣
e
i
θ
z=x+yi=\lvert z\rvert (\cos \theta + i\sin \theta)=\lvert z\rvert e^{i\theta}
z=x+yi=∣z∣(cosθ+isinθ)=∣z∣eiθ
单位根：
n
n
n次幂为
1
1
1的复数。
z
n
=
1
(
n
=
1
,
2
,
3
,
⋯
)
z^n=1 (n=1,2,3,\cdots)
zn=1(n=1,2,3,⋯)，这个方程的复数根
z
z
z为
n
n
n次单位根。
n
n
n次单位根有
n
n
n个，表达式为：
z
k
=
e
2
π
k
i
n
(
k
=
0
,
1
,
2
,
⋯
,
n
−
1
)
z_k=e^{\frac{2\pi ki}{n}}(k=0,1,2,\cdots,n-1)
zk=en2πki(k=0,1,2,⋯,n−1)
如果不加说明，一般叙述的
n
n
n次单位根，是指从
1
1
1开始逆时针方向的第一个解（不是
1
1
1），用
ω
n
\omega _n
ωn表示，其他解均用
ω
n
k
\omega_n^k
ωnk表示。
将单位根的概念放入图像中，则
w
n
w_n
wn与原点连线后，所得线与
x
x
x轴正半轴的夹角为
2
π
n
\frac{2\pi}{n}
n2π。
从图像中也可以看出，偶数次单位根有如下性质：
ω
n
i
=
−
ω
n
i
+
n
/
2
\omega_n^i=-\omega_n^{i+n/2}
ωni=−ωni+n/2。
2.快速傅里叶变换FFT
快速傅里叶变换(Fast Fourier Transform, FFT)，可以应用在两个
n
n
n度的多项式的乘法。相比于朴素的
O
(
n
2
)
O(n^2)
O(n2)做法，快速傅里叶变换有着更优秀的时间复杂度，是
O
(
n
log
n
)
O(n\log n)
O(nlogn)的。
由于两个整数的乘法也可以被当作多项式乘法，因此这个算法也可以用来加速大整数的乘法计算（令多项式中的
x
=
10
x=10
x=10）。
傅里叶变换
（简单了解即可）
傅里叶变换是一种分析信号的方法。对于一种连续性的信号，傅里叶变换可以分析信号的成分，也可以用这些成分合成信号。许多波形可以作为信号的成分。
设 f ( t ) f(t) f(t)是一个关于时间 t t t的函数，则傅里叶变换可以检测频率 ω \omega ω的周期在 f ( t ) f(t) f(t)的出现程度。
F ( ω ) = F [ f ( t ) ] = ∫ − ∞ ∞ f ( t ) e − i ω t d t F(\omega)=\mathbb{F}[f(t)]=\int_{-\infty}^{\infty} f(t)e^{-i\omega t} dt F(ω)=F[f(t)]=∫−∞∞f(t)e−iωtdt
逆变换为：
f ( t ) = F − 1 [ F ( ω ) ] = 1 2 π ∫ − ∞ ∞ F ( ω ) e i ω t d t f(t)=\mathbb{F}^{-1}[F(\omega)]=\frac{1}{2\pi }\int_{-\infty}^{\infty} F(\omega)e^{i\omega t} dt f(t)=F−1[F(ω)]=2π1∫−∞∞F(ω)eiωtdt
离散傅里叶变换DFT
先看一个数学问题：
计算 C n 3 + C n 7 + C n 11 + C n 15 + ⋯ C_n^3+C_n^7+C_n^{11}+C_n^{15}+\cdots Cn3+Cn7+Cn11+Cn15+⋯。
利用多项式求解这个问题，定义函数
f
(
x
)
f(x)
f(x)：
f
(
x
)
=
(
1
+
x
)
n
=
C
n
0
x
0
+
C
n
1
x
1
+
C
n
2
x
2
+
C
n
3
x
3
+
⋯
f(x)=(1+x)^n=C_n^0x^0+C_n^1x^1+C_n^2x^2+C_n^3x^3+\cdots
f(x)=(1+x)n=Cn0x0+Cn1x1+Cn2x2+Cn3x3+⋯
我们发现， f ( i ) = ( 1 + i ) n = C n 0 + C n 1 i − C n 2 − C n 3 i + ⋯ f(i)=(1+i)^n=C_n^0+C_n^1i-C_n^2-C_n^3i+\cdots f(i)=(1+i)n=Cn0+Cn1i−Cn2−Cn3i+⋯
尝试用这个式子凑出所求式，得到：
f
(
1
)
+
i
f
(
i
)
−
f
(
−
1
)
−
i
f
(
−
i
)
=
4
C
n
3
+
4
C
n
7
+
4
C
n
11
+
4
C
n
15
+
⋯
f(1)+if(i)-f(-1)-if(-i)=4C_n^3+4C_n^7+4C_n^{11}+4C_n^{15}+\cdots
f(1)+if(i)−f(−1)−if(−i)=4Cn3+4Cn7+4Cn11+4Cn15+⋯
所以，我们可得
C
n
3
+
C
n
7
+
C
n
11
+
C
n
15
+
⋯
=
2
n
+
i
(
1
+
i
)
n
−
i
(
1
−
i
)
n
4
C_n^3+C_n^7+C_n^{11}+C_n^{15}+\cdots=\frac{2^n+i(1+i)^n-i(1-i)^n}{4}
Cn3+Cn7+Cn11+Cn15+⋯=42n+i(1+i)n−i(1−i)n
解决这个数学问题的过程中，我们做了一件什么事情？
定义了一个连续函数 f ( x ) f(x) f(x)，然后赋予不同的 x x x，并通过和一些常值的四则运算，得到了我们要求的值。我们采用了 4 4 4次单位根作为 x x x值，拟合出了我们要求的答案。那么离散傅里叶变换的特点，是否也和单位根有关系呢？
傅里叶变换中，我们讨论的是连续的信号，在离散傅里叶变换中，我们要处理的信号就变成了一个个离散的值。如果说傅里叶变换是积分形式，那么离散傅里叶变换就是求和形式。
设
{
x
n
}
\{ x_n\}
{xn}是某一个满足有限性条件的序列，将DFT定义为：
X
k
=
∑
n
=
0
N
−
1
x
n
e
−
2
π
i
k
n
N
X_k=\displaystyle\sum\limits_{n=0}^{N-1} x_ne^{\frac{-2\pi ikn}{N}}
Xk=n=0∑N−1xneN−2πikn
这个形式和傅里叶变换中的积分形式比较相像。其逆变换表示为：
x n = ∑ n = 0 N − 1 X k e 2 π i k n N x_n=\displaystyle\sum\limits_{n=0}^{N-1} X_ke^{\frac{2\pi ikn}{N}} xn=n=0∑N−1XkeN2πikn
这里面，如果我们把 x n x_n xn看作多项式 f ( x ) f(x) f(x)的 x n x^n xn项的系数，这时计算得到的 X k X_k Xk恰好是 f ( x ) f(x) f(x)代入单位根 e − 2 π i k N e^{\frac{-2\pi ik}{N}} eN−2πik的点值 f ( e − 2 π i k N ) f(e^{\frac{-2\pi ik}{N}}) f(eN−2πik)。
所以离散傅里叶变换的特点，就是利用单位根 ω n \omega_n ωn解决问题，在多项式中的单位根处进行处理。
快速傅里叶变换FFT
FTT的基本思想是分治，就离散傅里叶变换DFT来说，FFT分治地求当 x = ω n k x=\omega_n^k x=ωnk时的 f ( x ) f(x) f(x)值。这个分治思想主要体现在将多项式分为奇次项和偶次项来处理。
原理
举例： f ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + a 4 x 4 + a 5 x 5 + a 6 x 6 + a 7 x 7 f(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7 f(x)=a0+a1x+a2x2+a3x3+a4x4+a5x5+a6x6+a7x7
按照次数的奇偶分成两组，变换成下面这种形式：
f
(
x
)
=
(
a
0
+
a
2
x
2
+
a
4
x
4
+
a
6
x
6
)
+
x
(
a
1
+
a
3
x
2
+
a
5
x
4
+
a
7
x
6
)
f(x)=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6)
f(x)=(a0+a2x2+a4x4+a6x6)+x(a1+a3x2+a5x4+a7x6)
分别用奇偶次项建立新的函数：
G
(
x
)
=
a
0
+
a
2
x
+
a
4
x
2
+
a
6
x
3
G(x)=a_0+a_2x+a_4x^2+a_6x^3
G(x)=a0+a2x+a4x2+a6x3
H
(
x
)
=
a
1
+
a
3
x
+
a
5
x
2
+
a
7
x
3
H(x)=a_1+a_3x+a_5x^2+a_7x^3
H(x)=a1+a3x+a5x2+a7x3
则
f
(
x
)
=
G
(
x
2
)
+
x
×
H
(
x
2
)
f(x)=G(x^2)+x\times H(x^2)
f(x)=G(x2)+x×H(x2)
偶数次单位根有如下性质： ω n i = − ω n i + n / 2 \omega_n^i=-\omega_n^{i+n/2} ωni=−ωni+n/2。同时我们注意到， G ( x ) G(x) G(x)和 H ( x ) H(x) H(x)都是偶函数（因为他们的项都是偶数次项），所以复平面上 ω n i \omega_n^i ωni和 ω n i + n / 2 \omega_n^{i+n/2} ωni+n/2的 G ( x 2 ) G(x^2) G(x2)和 H ( x 2 ) H(x^2) H(x2)对应的值应当是相同的。于是得到：
f ( ω n k ) = G ( ( ω n k ) 2 ) + ω n k × H ( ( ω n k ) 2 ) = G ( ω n 2 k ) + ω n k × H ( ω n 2 k ) f(\omega_n^k)=G((\omega_n^k)^2)+\omega_n^k\times H((\omega_n^k)^2)=G(\omega_n^{2k})+\omega_n^k\times H(\omega_n^{2k}) f(ωnk)=G((ωnk)2)+ωnk×H((ωnk)2)=G(ωn2k)+ωnk×H(ωn2k)
因为是偶数次根，从复平面上画图，不难得到 ω n 2 k = ω n / 2 k \omega _n^{2k}=\omega_{n/2}^k ωn2k=ωn/2k。所以继续推得：
f ( ω n k ) = G ( ω n / 2 k ) + ω n k × H ( ω n / 2 k ) f(\omega_n^k)=G(\omega_{n/2}^k)+\omega_n^k\times H(\omega_{n/2}^k) f(ωnk)=G(ωn/2k)+ωnk×H(ωn/2k)
同理得到：
f
(
ω
n
k
+
n
/
2
)
=
G
(
ω
n
2
k
+
n
)
+
ω
n
k
+
n
/
2
×
H
(
ω
n
2
k
+
n
)
f(\omega_n^{k+n/2})=G(\omega_n^{2k+n})+\omega_n^{k+n/2}\times H(\omega_n^{2k+n})
f(ωnk+n/2)=G(ωn2k+n)+ωnk+n/2×H(ωn2k+n)
利用
ω
n
i
=
−
ω
n
i
+
n
/
2
\omega_n^i=-\omega_n^{i+n/2}
ωni=−ωni+n/2变形，得到
f
(
ω
n
k
+
n
/
2
)
=
G
(
ω
n
2
k
+
n
)
−
ω
n
k
×
H
(
ω
n
2
k
+
n
)
f(\omega_n^{k+n/2})=G(\omega_n^{2k+n})-\omega_n^k\times H(\omega_n^{2k+n})
f(ωnk+n/2)=G(ωn2k+n)−ωnk×H(ωn2k+n)
在复平面上画图，不难得到
n
n
n次单位根以
n
n
n为周期，所以有
ω
n
k
=
ω
n
k
+
n
\omega_n^k=\omega_n^{k+n}
ωnk=ωnk+n，将其代入，得到
f
(
ω
n
k
+
n
/
2
)
=
G
(
ω
n
2
k
)
−
ω
n
k
×
H
(
ω
n
2
k
)
f(\omega_n^{k+n/2})=G(\omega_n^{2k})-\omega_n^k\times H(\omega_n^{2k})
f(ωnk+n/2)=G(ωn2k)−ωnk×H(ωn2k)
再利用
ω
n
2
k
=
ω
n
/
2
k
\omega _n^{2k}=\omega_{n/2}^k
ωn2k=ωn/2k，得到
f
(
ω
n
k
+
n
/
2
)
=
G
(
ω
n
/
2
k
)
−
ω
n
k
×
H
(
ω
n
/
2
k
)
f(\omega_n^{k+n/2})=G(\omega_{n/2}^k)-\omega_n^k\times H(\omega_{n/2}^k)
f(ωnk+n/2)=G(ωn/2k)−ωnk×H(ωn/2k)
现在的状况是，只要求出 G ( ω n / 2 k ) G(\omega_{n/2}^k) G(ωn/2k)和 H ( ω n / 2 k ) H(\omega_{n/2}^k) H(ωn/2k)，就可以直接求出 f ( ω n k ) f(\omega_n^k) f(ωnk)和 f ( ω n k + n / 2 ) f(\omega_n^{k+n/2}) f(ωnk+n/2)。
于是对 G G G和 H H H用同样的方法继续递归求解即可。
由于分治过程中处理的多项式长度是不断除以 2 2 2的，为了避免出现问题，我们通常需要让参与运算的多项式的长度满足 2 m 2^m 2m的形式，即多项式的度数是 2 2 2的幂数。如果长度不满足 2 2 2的幂数，需要将多项式在运算前补成长度为 2 m 2^m 2m的形式（补上的高次项的系数为 0 0 0即可），最后形成最高次项数为 2 m − 1 2^m-1 2m−1的多项式。
在代入值的时候，因为要代入 n n n个不同值，所以我们代入 ω n 0 , ω n 1 , ⋯ , ω n n − 1 \omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1} ωn0,ωn1,⋯,ωnn−1这 n = 2 m n=2^m n=2m个不同的值。
整个过程，将一个多项式从系数表示法变成了点值表示法。
位逆序变换 bit-reversal permutation
以 8 8 8项多项式为例，模拟拆分的过程如下：
在分治之前：
{
x
0
,
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
,
x
7
}
\{ x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7\}
{x0,x1,x2,x3,x4,x5,x6,x7}
第一次分治：
{
x
0
,
x
2
,
x
4
,
x
6
}
,
{
x
1
,
x
3
,
x
5
,
x
7
}
\{ x_0,x_2,x_4,x_6\},\{x_1,x_3,x_5,x_7\}
{x0,x2,x4,x6},{x1,x3,x5,x7}
第二次分治：
{
x
0
,
x
4
}
,
{
x
2
,
x
6
}
,
{
x
1
,
x
5
}
,
{
x
3
,
x
7
}
\{ x_0,x_4\},\{x_2,x_6\},\{x_1,x_5\},\{x_3,x_7\}
{x0,x4},{x2,x6},{x1,x5},{x3,x7}
第三次分治：
{
x
0
}
,
{
x
4
}
,
{
x
2
}
,
{
x
6
}
,
{
x
1
}
,
{
x
5
}
,
{
x
3
}
,
{
x
7
}
\{ x_0\},\{x_4\},\{x_2\},\{x_6\},\{x_1\},\{x_5\},\{x_3\},\{x_7\}
{x0},{x4},{x2},{x6},{x1},{x5},{x3},{x7}
将初始状态和最后的状态的二进制写出来看一看：
初始状态：
000
,
001
,
010
,
011
,
100
,
101
,
110
,
111
000,001,010,011,100,101,110,111
000,001,010,011,100,101,110,111
最后状态：
000
,
100
,
010
,
110
,
001
,
101
,
011
,
111
000,100,010,110,001,101,011,111
000,100,010,110,001,101,011,111
发现了什么？初始状态和最后状态的每一位的二进制都是相互颠倒的关系。也就是说，我们先模仿递归，把这些系数在原数组中先进行一下拆分，然后再倍增地合并这些算出来的结果。
采用逐位反转的方法，每个数的变换结果可以在
O
(
n
log
n
)
O(n\log n)
O(nlogn)的时间内求出。
（需要反转的位置有
n
n
n，每一个位置要调整二进制中的
log
n
\log n
logn位）
但是这个方法可以进一步优化：
设长度为 l e n = 2 k len=2^k len=2k，这里 k k k就是二进制的长度。设 r e v [ x ] rev[x] rev[x]表示 x x x的 k k k位二进制数（高位补 0 0 0）反转后得到的结果。已知 r e v [ 0 ] = 0 rev[0]=0 rev[0]=0。
接下来从小到大求
r
e
v
[
x
]
rev[x]
rev[x]。
当求
r
e
v
[
x
]
rev[x]
rev[x]时，
r
e
v
[
⌊
x
2
⌋
]
rev[\lfloor \frac{x}{2}\rfloor]
rev[⌊2x⌋]一定已经求解完毕。
⌊
x
2
⌋
\lfloor \frac{x}{2}\rfloor
⌊2x⌋可以通过
x
x
x二进制上右移一位得来，所以先将
x
x
x右移一位，然后将得到的结果翻转，再右移一位。此时得到的是
x
x
x除了二进制下从右向左的第一位以外其他位的反转的最终结果（因为右移一位时被挤掉了）。
剩下的一位就考虑
x
x
x的奇偶，如果
x
x
x原来是奇数，那么将得到的结果的最高位置
1
1
1，否则置
0
0
0。
举例来说，我们想翻转 ( 11001 ) 2 (11001)_2 (11001)2，首先第一步右移，得到 ( 01100 ) 2 (01100)_2 (01100)2，得到 ⌊ x 2 ⌋ \lfloor \frac{x}{2}\rfloor ⌊2x⌋。然后第二步翻转，得到 ( 00110 ) 2 (00110)_2 (00110)2。第三步再右移，得到 ( 00011 ) 2 (00011)_2 (00011)2。因为原来是奇数，所以最高位置 1 1 1，得到 ( 10011 ) 2 (10011)_2 (10011)2。这就是我们想要的翻转结果。
模拟分治（逆过程）
重新看一下我们分治的过程：
在分治之前：
{
x
0
,
x
1
,
x
2
,
x
3
,
x
4
,
x
5
,
x
6
,
x
7
}
\{ x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7\}
{x0,x1,x2,x3,x4,x5,x6,x7}
第一次分治：
{
x
0
,
x
2
,
x
4
,
x
6
}
,
{
x
1
,
x
3
,
x
5
,
x
7
}
\{ x_0,x_2,x_4,x_6\},\{x_1,x_3,x_5,x_7\}
{x0,x2,x4,x6},{x1,x3,x5,x7}
第二次分治：
{
x
0
,
x
4
}
,
{
x
2
,
x
6
}
,
{
x
1
,
x
5
}
,
{
x
3
,
x
7
}
\{ x_0,x_4\},\{x_2,x_6\},\{x_1,x_5\},\{x_3,x_7\}
{x0,x4},{x2,x6},{x1,x5},{x3,x7}
第三次分治：
{
x
0
}
,
{
x
4
}
,
{
x
2
}
,
{
x
6
}
,
{
x
1
}
,
{
x
5
}
,
{
x
3
}
,
{
x
7
}
\{ x_0\},\{x_4\},\{x_2\},\{x_6\},\{x_1\},\{x_5\},\{x_3\},\{x_7\}
{x0},{x4},{x2},{x6},{x1},{x5},{x3},{x7}
现在我们通过位逆序变换，把原来的顺序变成了第三次分治之后的顺序了。现在，我们就倒着从第三次分治的状态一步一步往回走。
逆过程就是：
以某次分治的逆过程为例，当前步骤是要将两个项数为
h
2
\frac{h}{2}
2h的多项式合并成一个长度为
h
h
h的多项式。
这时的单位根
ω
n
=
cos
(
2
π
h
)
+
i
sin
(
2
π
h
)
\omega_n=\cos(\frac{2\pi}{h})+i\sin(\frac{2\pi}{h})
ωn=cos(h2π)+isin(h2π)。
对于每一对相邻的项数为
h
2
\frac{h}{2}
2h的多项式，用下面两个公式将
f
(
ω
n
k
)
f(\omega_n^k)
f(ωnk)和
f
(
ω
n
k
+
n
/
2
)
f(\omega_n^{k+n/2})
f(ωnk+n/2)求解出来。（从左到右的
k
k
k依次为
0
,
1
,
2
,
⋯
0,1,2,\cdots
0,1,2,⋯）。
f
(
ω
n
k
)
=
G
(
ω
n
/
2
k
)
+
ω
n
k
×
H
(
ω
n
/
2
k
)
f(\omega_n^k)=G(\omega_{n/2}^k)+\omega_n^k\times H(\omega_{n/2}^k)
f(ωnk)=G(ωn/2k)+ωnk×H(ωn/2k)
f
(
ω
n
k
+
n
/
2
)
=
G
(
ω
n
/
2
k
)
−
ω
n
k
×
H
(
ω
n
/
2
k
)
f(\omega_n^{k+n/2})=G(\omega_{n/2}^k)-\omega_n^k\times H(\omega_{n/2}^k)
f(ωnk+n/2)=G(ωn/2k)−ωnk×H(ωn/2k)
快速傅里叶逆变换
我们用快速傅里叶变换，将多项式从系数表达式变成了点值表达式，从用一个个系数表示变成了单位根的形式来表示。但是我们最终要的形式还是一个系数的形式，所以需要再进行一次与快速傅里叶变换相逆的过程。
这种逆变换可以用正向的变换来表示。这里有两种理解的方式，但是线性代数角度比较好理解一点，所以这里只写线性代数的角度。
线性代数角度
离散傅里叶变换DFT本身是一个线性变换，我们可以理解为将目标多项式当作向量，左乘一个矩阵得到一个变换后的向量，以模拟把单位复根代入多项式的过程：
[ y 0 y 1 y 2 y 3 ⋮ y n − 1 ] = [ 1 1 1 1 ⋯ 1 1 ω n 1 ω n 2 ω n 3 ⋯ ω n n − 1 1 ω n 2 ω n 4 ω n 6 ⋯ ω n 2 ( n − 1 ) 1 ω n 3 ω n 6 ω n 9 ⋯ ω n 3 ( n − 1 ) ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ 1 ω n n − 1 ω n 2 ( n − 1 ) ω n 3 ( n − 1 ) ⋯ ω n ( n − 1 ) 2 ] [ a 0 a 1 a 2 a 3 ⋮ a n − 1 ] \left[ \begin{matrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots\\ y_{n-1} \end{matrix} \right]= \left[ \begin{matrix} 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)^2} \end{matrix} \right] \left[ \begin{matrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots\\ a_{n-1} \end{matrix} \right] ⎣ ⎡y0y1y2y3⋮yn−1⎦ ⎤=⎣ ⎡1111⋮11ωn1ωn2ωn3⋮ωnn−11ωn2ωn4ωn6⋮ωn2(n−1)1ωn3ωn6ωn9⋮ωn3(n−1)⋯⋯⋯⋯⋱⋯1ωnn−1ωn2(n−1)ωn3(n−1)⋮ωn(n−1)2⎦ ⎤⎣ ⎡a0a1a2a3⋮an−1⎦ ⎤
我们通过正向的快速傅里叶变换可以得到 y y y矩阵，这些 y y y就是点值表达式，即将各种单位根代入到 x x x当中得到的结果。
这个是正变换的矩阵乘法，那么逆变换，其实就是这个式子左右同时左乘中间的矩阵的逆阵即可。
然而这个矩阵有一个非常神奇的性质，那就是把每一项取倒数，然后再除以变换的长度，就可以得到这个矩阵的逆阵。（这里的长度指的是扩充长度后得到的 2 2 2的幂数的那个长度）
为了使计算的结果变为原来的倒数，我们可以得到：
1 ω k = ω k − 1 = e − 2 π i 2 = cos ( 2 π k ) + i sin ( − 2 π k ) \frac{1}{\omega_k}=\omega_k^{-1}=e^{-\frac{2\pi i}{2}}=\cos(\frac{2\pi}{k})+i \sin(-\frac{2\pi}{k}) ωk1=ωk−1=e−22πi=cos(k2π)+isin(−k2π)
注意到，我们只需要在正向的DFT中，把我们用的单位根 ω k \omega_k ωk变成 e − 2 π i k e^{-\frac{2\pi i}{k}} e−k2πi，计算结果就会变成原来的倒数了。唯一多出来的操作，就是将最后所有的结果除以过程中多项式的长度（是 2 2 2的幂数，不是多项式原本的项数）。其他的操作，和正向的DFT完全相同。
所以，我们可以定义一个函数，同时实现正向的快速傅里叶变换和反向的快速傅里叶逆变换，在参数中传入一个状态，表示采用正向变换还是采用逆向变换。
实现
1.复数
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x; y = _y;
}
Complex operator - (const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator + (const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator * (const Complex &b) const {
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
包含复数的加减乘运算。
2.位逆序变换
LL rev[4000005];
void change(Complex *y, int len) {
for (int i = 0; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1;
if (i & 1) rev[i] |= len >> 1;
}
for (int i = 0; i < len; ++i) {
if (i < rev[i]) swap(y[i], y[rev[i]]);
}
}
3.FFT
//on == 1: DFT; on == -1: IDFT
void FFT(Complex *y, int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j += h) {
Complex w(1, 0);
for (int k = j; k < j + h / 2; ++k) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1) {
for (int i = 0; i < len; ++i) {
y[i].x = y[i].x / len + 0.5;
}
}
}
LuoguP3803 - 【模板】多项式乘法（FFT）
先将系数读入，然后扩充成长度为 2 m 2^m 2m的形式。
步骤：
- 读入 n , m n,m n,m，根据 n , m n,m n,m确定最终的多项式项数 l e n len len。
- 读入第一个多项式，剩余项补0。
- 读入第二个多项式，剩余项补0。
- 将第一个多项式从系数表示法变为点值表示法（DFT）
- 将第二个多项式从系数表示法变为点值表示法（DFT）
- 将第二个多项式的点值乘到第一个多项式当中去。
- 将乘完后的第一个多项式从点值表示法变为系数表示法（IDFT）
- 输出实部作为答案
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const double PI = acos(-1.0);
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x; y = _y;
}
Complex operator - (const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator + (const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator * (const Complex &b) const {
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
LL rev[4000005];
void change(Complex *y, int len) {
for (int i = 0; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1;
if (i & 1) rev[i] |= len >> 1;
}
for (int i = 0; i < len; ++i) {
if (i < rev[i]) swap(y[i], y[rev[i]]);
}
}
void FFT(Complex *y, int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j += h) {
Complex w(1, 0);
for (int k = j; k < j + h / 2; ++k) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1) {
for (int i = 0; i < len; ++i) {
y[i].x = y[i].x / len + 0.5;
}
}
}
Complex x1[3000005], x2[3000005];
int n, m;
void main2() {
cin >> n >> m;
int len = 1;
while (len <= (n + m)) len <<= 1;
for (int i = 0; i <= n; ++i) {
cin >> x1[i].x;
x1[i].y = 0;
}
for (int i = n + 1; i < len; ++i) {
x1[i].x = x1[i].y = 0;
}
for (int i = 0; i <= m; ++i) {
cin >> x2[i].x;
x2[i].y = 0;
}
for (int i = m + 1; i < len; ++i) {
x2[i].x = x2[i].y = 0;
}
FFT(x1, len, 1); FFT(x2, len, 1);
for (int i = 0; i < len; ++i) {
x1[i] = x1[i] * x2[i];
}
FFT(x1, len, -1);
for (int i = 0; i <= n + m; ++i) {
cout << (LL)x1[i].x << ' ';
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
LL _ = 1;
// cin >> _;
while (_--) main2();
return 0;
}
“3变2”优化
回顾正常的多项式乘法，我们一共进行了 3 3 3次傅里叶变换，其中包括 2 2 2次正向变换， 1 1 1次逆向变换。
但是每一次都是 O ( n log n ) O(n\log n) O(nlogn)，进行 3 3 3次，如果 n n n稍微大一点，这个常数就会觉得稍微有点大了。
这里我们有一个优化，先看下面这个式子：
(
a
+
b
i
)
2
=
a
2
−
b
2
+
2
a
b
i
(a+bi)^2=a^2-b^2+2abi
(a+bi)2=a2−b2+2abi
于是，我们可以进行 f ( x ) ⋅ g ( x ) f(x)\cdot g(x) f(x)⋅g(x)时，直接把 g ( x ) g(x) g(x)放到 f ( x ) f(x) f(x)的虚部，然后求 f ( x ) 2 f(x)^2 f(x)2，然后把 f ( x ) f(x) f(x)的虚部取出来除以 2 2 2就是答案了。
LuoguP3803 - 【模板】多项式乘法（FFT）
“3变2”优化代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const double PI = acos(-1.0);
struct Complex {
double x, y;
Complex(double _x = 0.0, double _y = 0.0) {
x = _x; y = _y;
}
Complex operator - (const Complex &b) const {
return Complex(x - b.x, y - b.y);
}
Complex operator + (const Complex &b) const {
return Complex(x + b.x, y + b.y);
}
Complex operator * (const Complex &b) const {
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
LL rev[4000005];
void change(Complex *y, int len) {
for (int i = 0; i < len; ++i) {
rev[i] = rev[i >> 1] >> 1;
if (i & 1) rev[i] |= len >> 1;
}
for (int i = 0; i < len; ++i) {
if (i < rev[i]) swap(y[i], y[rev[i]]);
}
}
void FFT(Complex *y, int len, int on) {
change(y, len);
for (int h = 2; h <= len; h <<= 1) {
Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));
for (int j = 0; j < len; j += h) {
Complex w(1, 0);
for (int k = j; k < j + h / 2; ++k) {
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1) {
for (int i = 0; i < len; ++i) {
y[i].x = y[i].x / len + 0.5;
y[i].y = y[i].y / (len * 2) + 0.5;
}
}
}
Complex x1[3000005];
int n, m;
void main2() {
cin >> n >> m;
int len = 1;
while (len <= (n + m)) {
len <<= 1;
}
for (int i = 0; i <= n; ++i) {
cin >> x1[i].x;
}
for (int i = n + 1; i < len; ++i) {
x1[i].x = 0;
}
for (int i = 0; i <= m; ++i) {
cin >> x1[i].y;
}
for (int i = m + 1; i < len; ++i) {
x1[i].y = 0;
}
FFT(x1, len, 1);
for (int i = 0; i < len; ++i) {
x1[i] = x1[i] * x1[i];
}
FFT(x1, len, -1);
for (int i = 0; i <= n + m; ++i) {
cout << (LL)x1[i].y << ' ';
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0); cout.tie(0);
LL _ = 1;
// cin >> _;
while (_--) main2();
return 0;
}
1958
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
