O ( n log n ) O(n\log n) O(nlogn)的时间将一个多项式的系数表达式变成点值表达式（这两种表达方式后文会提到）。
O ( n log n ) O(n\log n) O(nlogn)的时间将一个多项式的点值表达式变成系数表达式（这两种表达方式后文会提到）。
O ( n log n ) O(n\log n) O(nlogn)的时间求两个多项式的乘积。
一个
n
n
n项的整式多项式一定可以表示成：
f
(
x
)
=
∑
i
=
0
n
−
1
(
a
i
x
i
)
f(x)=\sum\limits_{i=0}^{n-1}(a_ix^i)
f(x)=i=0∑n−1(aixi)
（有些
a
i
a_i
ai可能为
0
0
0）
所以，一个系数
a
a
a的（有序）集合可以唯一确定一个多项式。
例如：
a
=
{
1
,
2
,
0
,
3
,
5
}
a=\{1,2,0,3,5\}
a={1,2,0,3,5}表示的就是多项式：
f
(
x
)
=
1
+
2
x
+
3
x
3
+
5
x
4
f(x)=1+2x+3x^3+5x^4
f(x)=1+2x+3x3+5x4。
这里
a
=
{
1
,
2
,
0
,
3
,
5
}
a=\{1,2,0,3,5\}
a={1,2,0,3,5}叫
f
(
x
)
f(x)
f(x)的系数表达式。
所以，所谓求两个多项式的乘积，就是已知两个多项式的系数表达式，求它们乘积的多项式的系数表达式。
（既然大家都这样说，我也这样说吧）
之前说的 f ( x ) = ∑ i = 0 n − 1 ( a i x i ) f(x)=\sum\limits_{i=0}^{n-1}(a_ix^i) f(x)=i=0∑n−1(aixi)中的 n n n就是它的阶，注意 n n n阶多项式不一定是 n n n次的，因为 a n − 1 a_{n-1} an−1可能为 0 0 0。
前面已说。
一个 n n n次多项式就是一个函数，取它图像上的 n + 1 n+1 n+1个不同的点就可以确定这个多项式。
例如 f ( x ) f(x) f(x)是个 4 4 4次多项式，取它的图像上的任意 5 5 5个不同的点： ( x 1 , f ( x 1 ) ) (x_1,f(x_1)) (x1,f(x1))， ( x 2 , f ( x 2 ) ) (x_2,f(x_2)) (x2,f(x2))，…， ( x 5 , f ( x 5 ) ) (x_5,f(x_5)) (x5,f(x5))，可以确定这个多项式的系数，因为可以列出一个 5 5 5元一次方程组解出系数。
于是一个
n
n
n次多项式的点值表达式就是它的图像上
n
+
1
n+1
n+1个不同的点。
换句话说，
n
+
1
n+1
n+1个不同的点可以唯一确定一个
n
n
n次多项式。
（以上两句话式FFT和逆FFT的核心）
在实数范围内，老师会告诉你负数没有平方根，所以我们为了让负数有平方根，将数的范围从实数扩展到了复数。
（实数也属于复数，复数另外一部分是虚数）
引入一个符号
i
i
i，
i
2
=
−
1
i^2=-1
i2=−1。
于是任意一个复数都可以表示成
a
+
b
i
a+bi
a+bi，其中
a
a
a和
b
b
b都是实数。
例如：
x
2
=
−
7
x^2=-7
x2=−7的根是：
x
1
=
7
i
x_1=\sqrt 7i
x1=7
i，
x
2
=
−
7
i
x_2=-\sqrt 7i
x2=−7
i。
加减乘都和整式的运算是一样的，例如两个复数
z
1
=
a
1
+
b
1
i
z_1=a_1+b_1i
z1=a1+b1i，
z
2
=
a
2
+
b
2
i
z_2=a_2+b_2i
z2=a2+b2i相乘：
z
1
z
2
=
(
a
1
+
b
1
i
)
(
a
2
+
b
2
i
)
=
(
a
1
a
2
−
b
1
b
2
)
+
(
a
1
b
2
+
a
2
b
1
)
i
z_1z_2=(a_1+b_1i)(a_2+b_2i)=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i
z1z2=(a1+b1i)(a2+b2i)=(a1a2−b1b2)+(a1b2+a2b1)i
除法这里不会用到，但也不难，可以分母实数化：
来自百度百科：复数除法
数轴可以表示一切实数，平面可以表示一切复数：
（其实这个图没有什么用）
在复平面内的一个点
(
x
,
y
)
(x,y)
(x,y)表示一个复数
x
+
y
i
x+yi
x+yi。
定义：
ω
n
k
=
c
o
s
(
2
π
n
k
)
+
s
i
n
(
2
π
n
k
)
i
\omega_n^k=cos\left(\frac{2\pi}{n}k\right)+sin\left(\frac{2\pi}{n}k\right)i
ωnk=cos(n2πk)+sin(n2πk)i
ω
n
1
\omega_n^1
ωn1叫
n
n
n次单位复根。
很难懂是不是？把
ω
8
0
\omega_8^0
ω80，…，
ω
8
7
\omega_8^7
ω87对应的点在复平面上画出来你就知道了。
（由于几何画板的公式编辑太弱，所以只标了
A
B
C
ABC
ABC）
点
A
A
A，
B
B
B，
C
C
C…对应的复数分别是
ω
8
0
\omega_8^0
ω80，
ω
8
1
\omega_8^1
ω81，
ω
8
2
\omega_8^2
ω82，…
为什么？数学必修四的单位圆与三角函数会告诉你的。
证明：
由于 e i θ = c o s θ + s i n θ ⋅ i e^{i\theta}=cos\theta+sin\theta\cdot i eiθ=cosθ+sinθ⋅i（我不会证，貌似要用泰勒展开）
所以 ω n i = e 2 π n i \omega_n^i=e^{\frac{2\pi}{n}i} ωni=en2πi
所以 ω n i ω n j = e 2 π n i e 2 π n j = e 2 π n ( i + j ) = ω n i + j \omega_n^i\omega_n^j=e^{\frac{2\pi}{n}i}e^{\frac{2\pi}{n}j}=e^{\frac{2\pi}{n}(i+j)}=\omega_n^{i+j} ωniωnj=en2πien2πj=en2π(i+j)=ωni+j
证明：
ω 2 n 2 k = c o s ( 2 π 2 n 2 k ) + s i n ( 2 π 2 n 2 k ) i = c o s ( 2 π n k ) + s i n ( 2 π n k ) i = ω n k \omega_{2n}^{2k}=cos\left(\frac{2\pi}{2n}2k\right)+sin\left(\frac{2\pi}{2n}2k\right)i=cos\left(\frac{2\pi}{n}k\right)+sin\left(\frac{2\pi}{n}k\right)i=\omega_n^k ω2n2k=cos(2n2π2k)+sin(2n2π2k)i=cos(n2πk)+sin(n2πk)i=ωnk
证明： ω n k + n 2 = c o s ( 2 π n k + π ) + s i n s ( 2 π n k + π ) i = − c o s ( 2 π n k ) − s i n ( 2 π n k ) i = − ω n k \omega_{n}^{k+\frac{n}{2}}=cos\left(\frac{2\pi}{n}k+\pi\right)+sins\left(\frac{2\pi}{n}k+\pi\right)i=-cos\left(\frac{2\pi}{n}k\right)-sin\left(\frac{2\pi}{n}k\right)i=-\omega_n^k ωnk+2n=cos(n2πk+π)+sins(n2πk+π)i=−cos(n2πk)−sin(n2πk)i=−ωnk
证明：同上一个证明
证明：
∑ i = 0 n − 1 ( ω n k ) i = 1 + ω n k + ( ω n k ) 2 + . . . + ( ω n k ) n − 1 = ( ω n k ) n − 1 ω n k − 1 = ω n k n − 1 ω n k − 1 = 1 − 1 ω n k − 1 = 0 \sum\limits_{i=0}^{n-1}(\omega_n^k)^i=1+\omega_n^k+(\omega_n^k)^2+...+(\omega_n^k)^{n-1}=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{\omega_n^{kn}-1}{\omega_n^k-1}=\frac{1-1}{\omega_n^k-1}=0 i=0∑n−1(ωnk)i=1+ωnk+(ωnk)2+...+(ωnk)n−1=ωnk−1(ωnk)n−1=ωnk−1ωnkn−1=ωnk−11−1=0
已知
k
1
k_1
k1阶多项式
f
(
x
)
f(x)
f(x)和
k
2
k_2
k2阶多项式
g
(
x
)
g(x)
g(x)的系数表达式。
我们要求的是
k
1
+
k
2
k_1+k_2
k1+k2阶多项式
h
(
x
)
=
f
(
x
)
⋅
g
(
x
)
h(x)=f(x)\cdot g(x)
h(x)=f(x)⋅g(x)的系数表达式。
忽略极大的常数，暴力的时间复杂度是
O
(
n
2
)
O(n^2)
O(n2)，FFT可以把它优化到
O
(
n
log
n
)
O(n\log n)
O(nlogn)。
FFT（IFFT）是通过计算一半的点值表达式得到另一半，要计算的那一半重复前面的步骤。
[特别提醒]
我们代入求点值表达式的 n n n个 x x x是： ω n 0 , . . . , ω n n − 1 \omega_n^0,...,\omega_n^{n-1} ωn0,...,ωnn−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+...+a_{n-1}x^{n-1}
A(x)=a0+a1x+a2x2+...+an−1xn−1
令
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+...+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+...+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)（可以自己带进去算一下）
那么
A
(
ω
n
k
)
=
A
0
(
ω
n
2
k
)
+
ω
n
k
A
1
(
ω
n
2
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_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k})
A(ωnk)=A0(ωn2k)+ωnkA1(ωn2k)=A0(ω2nk)+ωnkA1(ω2nk)
A
(
ω
n
k
+
n
2
)
=
A
0
(
ω
n
2
k
+
n
)
+
ω
n
k
+
n
2
A
1
(
ω
n
2
k
+
n
)
=
A
0
(
ω
n
2
k
)
−
ω
n
k
A
1
(
ω
n
2
k
)
=
A
0
(
ω
n
2
k
)
−
ω
n
k
A
1
(
ω
n
2
k
)
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_0(\omega_{n}^{2k})-\omega_n^kA_1(\omega_{n}^{2k})=A_0(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_1(\omega_{\frac{n}{2}}^{k})
A(ωnk+2n)=A0(ωn2k+n)+ωnk+2nA1(ωn2k+n)=A0(ωn2k)−ωnkA1(ωn2k)=A0(ω2nk)−ωnkA1(ω2nk)
这两个只有正负号不同！
我们定义FFT(A,n)，表示将系数表达式的系数集合
A
A
A转换为点值表达式的
y
y
y值集合（代入的
x
x
x是
ω
n
0
,
ω
n
1
,
.
.
.
ω
n
n
−
1
\omega_n^0,\omega_n^1,...\omega_n^{n-1}
ωn0,ωn1,...ωnn−1）。
那么，FFT(A,n)中，先提出
A
0
A_0
A0和
A
1
A_1
A1，再执行FFT(A0,n/2)和FFT(A1,n/2)。
执行完后的
A
0
A_0
A0和
A
1
A_1
A1已经变成了当 “ 代入的
x
x
x是
ω
n
2
0
,
ω
n
2
1
,
.
.
.
ω
n
2
n
2
−
1
\omega_\frac{n}{2}^0,\omega_\frac{n}{2}^1,...\omega_\frac{n}{2}^{\frac{n}{2}-1}
ω2n0,ω2n1,...ω2n2n−1 ” 时，对应的
y
y
y的集合。
然后就可以循环从
0
0
0到
n
2
\frac{n}{2}
2n，A[i]=A0[i]+w*A1[i]，A[i+n/2]=A0[i]-w*A1[i]。
等讲了IFFT的递归版一起看。
（来自不知名的大佬）
这是
n
=
8
n=8
n=8时的递归FFT的系数示意图。
发现最后得到的系数序列，每个数的二进制反转后恰好是最开始的序列（为什么？用心去感受一下吧）。
所以只需要把最开始的系数变成最后的样子，再一层层向上合并即可。
也就是说，把最开始的每个系数和它二进制反转后的数交换位置。
定义R[i]表示i的二进制反转后的数，L是i的二进制位数（即
log
2
n
\log_2 n
log2n），那么有：
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1))
下面是解释：
R[i>>1]表示i除去最后一位后前面几位的反转，那么把它和i的最后一位换个位置就好了。
等讲了IFFT的非递归版一起看。
( y 0 y 1 y 2 y 3 ⋮ y n − 1 ) = ( ω n 0 ω n 0 ω n 0 ω n 0 ⋯ ω n 0 ω n 0 ω n 1 ω n 2 ω n 3 ⋯ ω n n − 1 ω n 0 ω n 2 ω n 4 ω n 6 ⋯ ω n 2 ( n − 1 ) ω n 0 ω n 3 ω n 6 ω n 9 ⋯ ω n 3 ( n − 1 ) ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ω n 0 ω n n − 1 ω n 2 ( n − 1 ) ω n 3 ( n − 1 ) ⋯ ω n ( n − 1 ) ( n − 1 ) ) × ( 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} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1}\\ \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)}\\ \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}\\ \end{matrix} \right)\times\\ \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right) ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛ωn0ωn0ωn0ωn0⋮ωn0ωn0ωn1ωn2ωn3⋮ωnn−1ωn0ωn2ωn4ωn6⋮ωn2(n−1)ωn0ωn3ωn6ωn9⋮ωn3(n−1)⋯⋯⋯⋯⋱⋯ωn0ωnn−1ωn2(n−1)ωn3(n−1)⋮ωn(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞×⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
FFT可以表示为这样，
{
y
}
\{y\}
{y}是点值表达式的
y
y
y值集合，
{
a
}
\{a\}
{a}是系数表达式的系数集合。
（矩阵乘法可以自己百度一下，在这里就是
y
k
=
∑
i
=
0
n
−
1
(
a
i
ω
n
k
i
)
y_k=\sum\limits_{i=0}^{n-1}\left(a_i\omega_n^{ki}\right)
yk=i=0∑n−1(aiωnki)）
例如： y 2 = f ( x 2 ) = f ( ω n 2 ) = a 0 ( ω n 2 ) 0 + a 1 ( ω n 2 ) 1 + a 2 ( ω n 2 ) 2 + . . . + a n − 1 ( ω n 2 ) n − 1 y_2=f(x_2)=f(\omega_n^2)=a_0(\omega_n^2)^0+a_1(\omega_n^2)^1+a_2(\omega_n^2)^2+...+a_{n-1}(\omega_n^2)^{n-1} y2=f(x2)=f(ωn2)=a0(ωn2)0+a1(ωn2)1+a2(ωn2)2+...+an−1(ωn2)n−1 = a 0 ω n 0 + a 1 ω n 2 + a 2 ω n 4 + . . . + a n − 1 ω n 2 ( n − 1 ) =a_0\omega_n^0+a_1\omega_n^2+a_2\omega_n^4+...+a_{n-1}\omega_n^{2(n-1)} =a0ωn0+a1ωn2+a2ωn4+...+an−1ωn2(n−1) = ∑ i = 0 n − 1 ( a i ω n 2 i ) =\sum\limits_{i=0}^{n-1}\left(a_i\omega_n^{2i}\right) =i=0∑n−1(aiωn2i)
令
V
=
(
ω
n
0
ω
n
0
ω
n
0
ω
n
0
⋯
ω
n
0
ω
n
0
ω
n
1
ω
n
2
ω
n
3
⋯
ω
n
n
−
1
ω
n
0
ω
n
2
ω
n
4
ω
n
6
⋯
ω
n
2
(
n
−
1
)
ω
n
0
ω
n
3
ω
n
6
ω
n
9
⋯
ω
n
3
(
n
−
1
)
⋮
⋮
⋮
⋮
⋱
⋮
ω
n
0
ω
n
n
−
1
ω
n
2
(
n
−
1
)
ω
n
3
(
n
−
1
)
⋯
ω
n
(
n
−
1
)
(
n
−
1
)
)
V=\left( \begin{matrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1}\\ \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)}\\ \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}\\ \end{matrix} \right)
V=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛ωn0ωn0ωn0ωn0⋮ωn0ωn0ωn1ωn2ωn3⋮ωnn−1ωn0ωn2ωn4ωn6⋮ωn2(n−1)ωn0ωn3ωn6ωn9⋮ωn3(n−1)⋯⋯⋯⋯⋱⋯ωn0ωnn−1ωn2(n−1)ωn3(n−1)⋮ωn(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
则
(
y
0
y
1
y
2
y
3
⋮
y
n
−
1
)
=
V
×
(
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)= V\times\\ \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right)
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=V×⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
那么
(
y
0
y
1
y
2
y
3
⋮
y
n
−
1
)
×
V
−
1
=
(
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)\times V^{-1}= \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right)
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛y0y1y2y3⋮yn−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞×V−1=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a0a1a2a3⋮an−1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
也就是说，我们只要构造出 V − 1 V^{-1} V−1，IFFT就解决了。
换句话说，要构造一个矩阵
V
−
1
V^{-1}
V−1，使得
V
×
V
−
1
=
A
V\times V^{-1}=A
V×V−1=A（
A
A
A是单位矩阵）
其中
A
=
(
1
0
0
⋯
0
0
1
0
⋯
0
0
0
1
⋯
0
⋮
⋮
⋮
⋱
⋮
0
0
0
⋯
1
)
A= \left( \begin{matrix} 1 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1 \end{matrix} \right)
A=⎝⎜⎜⎜⎜⎜⎛100⋮0010⋮0001⋮0⋯⋯⋯⋱⋯000⋮1⎠⎟⎟⎟⎟⎟⎞
（计算一下会发现任何矩阵乘 A A A都是它本身）
于是傅里叶说，
V
−
1
=
(
1
n
ω
n
0
1
n
ω
n
0
1
n
ω
n
0
1
n
ω
n
0
⋯
1
n
ω
n
0
1
n
ω
n
0
1
n
ω
n
−
1
1
n
ω
n
−
2
1
n
ω
n
−
3
⋯
1
n
ω
n
−
(
n
−
1
)
1
n
ω
n
0
1
n
ω
n
−
2
1
n
ω
n
−
4
1
n
ω
n
−
6
⋯
1
n
ω
n
−
2
(
n
−
1
)
1
n
ω
n
0
1
n
ω
n
−
3
1
n
ω
n
−
6
1
n
ω
n
−
9
⋯
1
n
ω
n
−
3
(
n
−
1
)
⋮
⋮
⋮
⋮
⋱
⋮
1
n
ω
n
0
1
n
ω
n
−
(
n
−
1
)
1
n
ω
n
−
2
(
n
−
1
)
1
n
ω
n
−
3
(
n
−
1
)
⋯
1
n
ω
n
−
(
n
−
1
)
(
n
−
1
)
)
V^{-1}= \left( \begin{matrix} \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \cdots & \frac{1}{n}\omega_n^0\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-1} & \frac{1}{n}\omega_n^{-2} & \frac{1}{n}\omega_n^{-3} & \cdots & \frac{1}{n}\omega_n^{-(n-1)}\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-2} & \frac{1}{n}\omega_n^{-4} & \frac{1}{n}\omega_n^{-6} & \cdots & \frac{1}{n}\omega_n^{-2(n-1)}\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-3} & \frac{1}{n}\omega_n^{-6} & \frac{1}{n}\omega_n^{-9} & \cdots & \frac{1}{n}\omega_n^{-3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-(n-1)} & \frac{1}{n}\omega_n^{-2(n-1)} & \frac{1}{n}\omega_n^{-3(n-1)} & \cdots & \frac{1}{n}\omega_n^{-(n-1)(n-1)}\\ \end{matrix} \right)
V−1=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛n1ωn0n1ωn0n1ωn0n1ωn0⋮n1ωn0n1ωn0n1ωn−1n1ωn−2n1ωn−3⋮n1ωn−(n−1)n1ωn0n1ωn−2n1ωn−4n1ωn−6⋮n1ωn−2(n−1)n1ωn0n1ωn−3n1ωn−6n1ωn−9⋮n1ωn−3(n−1)⋯⋯⋯⋯⋱⋯n1ωn0n1ωn−(n−1)n1ωn−2(n−1)n1ωn−3(n−1)⋮n1ωn−(n−1)(n−1)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
简单点说，（下标从零开始） V i , j = ω n i j V_{i,j}=\omega_n^{ij} Vi,j=ωnij， V − 1 i , j = 1 n ω n − i j V^{-1}{}_{i,j}=\frac{1}{n}\omega_n^{-ij} V−1i,j=n1ωn−ij。
接下来要证明 V × V − 1 = A = ( 1 0 0 ⋯ 0 0 1 0 ⋯ 0 0 0 1 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ 1 ) V\times V^{-1}=A=\left( \begin{matrix} 1 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1 \end{matrix} \right) V×V−1=A=⎝⎜⎜⎜⎜⎜⎛100⋮0010⋮0001⋮0⋯⋯⋯⋱⋯000⋮1⎠⎟⎟⎟⎟⎟⎞
根据矩阵乘法的定义， A x , y = ∑ i = 0 n − 1 ( V x , i V − 1 i , y ) A_{x,y}=\sum\limits_{i=0}^{n-1}(V_{x,i}V^{-1}{}_{i,y}) Ax,y=i=0∑n−1(Vx,iV−1i,y) = ∑ i = 0 n − 1 ( ω n x i 1 n ω n − i y ) =\sum\limits_{i=0}^{n-1}\left(\omega_n^{xi}\frac{1}{n}\omega_n^{-iy}\right) =i=0∑n−1(ωnxin1ωn−iy) = 1 n ∑ i = 0 n − 1 ( ω n x − y ) i =\frac{1}{n}\sum\limits_{i=0}^{n-1}(\omega_n^{x-y})^i =n1i=0∑n−1(ωnx−y)i
当
x
=
y
x=y
x=y时，
ω
n
x
−
y
=
1
\omega_n^{x-y}=1
ωnx−y=1，
A
i
,
j
=
1
A_{i,j}=1
Ai,j=1；
当
x
≠
y
x\neq y
x=y时，根据求和引理，
A
i
,
j
=
0
A_{i,j}=0
Ai,j=0。
得证。
综上所述：
既然构造出来了
V
−
1
V^{-1}
V−1，我们就知道要代入计算的
x
x
x集合就是：
{
ω
n
0
,
ω
n
−
1
,
.
.
.
,
ω
n
−
(
n
−
1
)
}
\{\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}\}
{ωn0,ωn−1,...,ωn−(n−1)}
为了减少误差，最后统一乘 1 n \frac{1}{n} n1。
把 “ 求多项式乘积的基本步骤 ” 中的第3步得到的 h ( x ) h(x) h(x)的点值表达式的 y y y集合当做一个多项式的系数表达式集合，用一次FFT求出 x x x集合是 { ω n 0 , ω n − 1 , . . . , ω n − ( n − 1 ) } \{\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}\} {ωn0,ωn−1,...,ωn−(n−1)}时，这个多项式的点值表达式的 y y y集合，就是 h ( x ) h(x) h(x)的系数集合的 n n n倍。
还记得FFT算法流程的这一步吗：
A[i]=A0[i]+w*A1[i]，A[i+n/2]=A0[i]-w*A1[i]
将它变为这样就是IFFT了：
A[i]=A0[i]-w*A1[i]，A[i+n/2]=A0[i]+w*A1[i]
现在知道为什么代码要一起展示了，因为只需在FFT函数中加一个参数就可以实现逆变换了。
看下面。
大整数乘法一般的时间复杂度是 O ( n 2 ) O(n^2) O(n2)，把两个大整数看成两个多项式的系数，求它们的乘积，最后处理进位，就可以做到 O ( n log n ) O(n\log n) O(nlogn)了。
据说递归版要爆空间（我这道题和洛谷上的那道都用递归的过了），所以大家都写非递归版。
#include<cmath>
#include<cstdio>
#include<cstring>
#define LOG 17
#define MAXN (1<<LOG)
int res[MAXN+5];
char x[MAXN+5],y[MAXN+5];
struct Complex{
double x,y;//x+yi
Complex(){x=y=0;}
Complex(double a,double b):x(a),y(b){}
Complex operator + (Complex b)const{return Complex(x+b.x,y+b.y);}
Complex operator - (Complex b)const{return Complex(x-b.x,y-b.y);}
Complex operator * (Complex b)const{return Complex(x*b.x-y*b.y,y*b.x+x*b.y);}
}A0[MAXN+5],B0[MAXN+5];
//复数类，好像C++有一个叫complex的库
#define PI acos(-1)
void FFT(Complex *A,int len,int sign){
if(len==1)
return;
int mid=len>>1;
Complex A1[mid+5],A2[mid+5];//按奇偶分离系数到A1和A2
for(int i=0;i<mid;i++)
A1[i]=A[i<<1],A2[i]=A[i<<1|1];
FFT(A1,mid,sign),FFT(A2,mid,sign);
for(int i=0;i<mid;i++){
Complex w(cos(2.0*PI/len*i),sign*sin(2.0*PI/len*i));
//w是len次单位复根的i次方
//sign=-1是逆变换
A[i]=A1[i]+w*A2[i],A[i+mid]=A1[i]-w*A2[i];
}
}
int main(){
while(~scanf("%s%s",x+1,y+1)){
int len1=strlen(x+1),len2=strlen(y+1),N=1;
while(N<len1+len2)
N<<=1;
//统一长度为2的幂
for(int i=0;i<len1;i++)
A0[i]=Complex(x[len1-i]-'0',0);
for(int i=0;i<len2;i++)
B0[i]=Complex(y[len2-i]-'0',0);
//倒着存可以避免很多恶心的问题
FFT(A0,N,1);
FFT(B0,N,1);
for(int i=0;i<N;i++)
A0[i]=A0[i]*B0[i];
//得到乘积多项式的点值表达式
FFT(A0,N,-1);
for(int i=0;i<N;i++)
res[i+1]=int(A0[i].x/N+0.5);
//注意最后统一除以n，为避免精度误差，+0.5再取整
for(int i=1;i<=N;i++)
res[i+1]+=res[i]/10,res[i]%=10;
while(N>1&&res[N]==0)
N--;
//处理进位和前导零
for(int i=N;i>=1;i--)
putchar(res[i]+'0');
putchar('\n');
memset(A0,0,sizeof A0);
memset(B0,0,sizeof B0);
memset(res,0,sizeof res);
}
}
建议写这种。
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define LOG 17
#define MAXN (1<<LOG)
int res[MAXN+5];
char x[MAXN+5],y[MAXN+5];
struct Complex{
double x,y;//x+yi
Complex(){x=y=0;}
Complex(double a,double b):x(a),y(b){}
Complex operator + (Complex b)const{return Complex(x+b.x,y+b.y);}
Complex operator - (Complex b)const{return Complex(x-b.x,y-b.y);}
Complex operator * (Complex b)const{return Complex(x*b.x-y*b.y,y*b.x+x*b.y);}
}A0[MAXN+5],B0[MAXN+5];
#define PI acos(-1)
int R[MAXN+5];
void FFT(Complex *A,int len,int sign){
for(int i=0;i<len;i++)
if(i<R[i])//不加的话会swap两次，最后什么都没有发生
swap(A[i],A[R[i]]);
//调整成最后的顺序
for(int mid=1;mid<len;mid<<=1)//枚举当前要处理的区间长度的一半
for(int i=0;i<len;i+=mid<<1)//每2*mid的长度为一段，分开处理
for(int j=0;j<mid;j++){//处理当前2*mid长度的区间
Complex w(cos(PI/mid*j),sign*sin(PI/mid*j));//2.0被mid消掉了
Complex tmp1=A[i+j],tmp2=A[i+j+mid];
A[i+j]=tmp1+w*tmp2;
A[i+j+mid]=tmp1-w*tmp2;
}
}
int main(){
while(~scanf("%s%s",x+1,y+1)){
int len1=strlen(x+1),len2=strlen(y+1),N=1,Log=0;
while(N<len1+len2)
N<<=1,Log++;
for(int i=0;i<len1;i++)
A0[i]=Complex(x[len1-i]-'0',0);
for(int i=0;i<len2;i++)
B0[i]=Complex(y[len2-i]-'0',0);
for(int i=0;i<N;i++)
R[i]=(R[i>>1]>>1)|((i&1)<<(Log-1));
//dp得到R数组
FFT(A0,N,1),FFT(B0,N,1);
for(int i=0;i<N;i++)
A0[i]=A0[i]*B0[i];
FFT(A0,N,-1);
for(int i=0;i<N;i++)
res[i+1]=int(A0[i].x/N+0.5);
for(int i=1;i<=N;i++)
res[i+1]+=res[i]/10,res[i]%=10;
while(N>1&&res[N]==0)
N--;
for(int i=N;i>=1;i--)
putchar(res[i]+'0');
putchar('\n');
memset(A0,0,sizeof A0);
memset(B0,0,sizeof B0);
memset(res,0,sizeof res);
}
}
感谢这些大佬：