多项式算法——快速傅里叶变换FFT
目录
两个多项式函数相加的时间复杂度为 O ( n ) O(n) O(n),而相乘的时间复杂度为 O ( n 2 ) O(n^2) O(n2),即所有系数的笛卡尔积。傅里叶变换DFT将优化多项式相乘(以下简称卷积),而快速傅里叶变换则优化至 O ( n log n ) O(n \log n) O(nlogn)。
前置知识:
- 微积分(欧拉公式)
- 复变函数(单位元根)
- 线性代数(矩阵方程)
多项式
本章定义的多项式函数为定义在实数域上,关于变量 x x x形如:
A ( x ) = ∑ j = 0 n − 1 a j x j A(x) = \sum_{j=0}^{n-1} a_j x^j A(x)=j=0∑n−1ajxj
的 n n n次多项式, a j a_j aj为多项式的系数,如果多项式非零最高次数项的系数为 a k a_k ak,那么我们说 A ( x ) A(x) A(x)的次数为 k k k,记作 d e g r e e ( A ) = k degree(A)=k degree(A)=k。
任何一个严格大于多项式 A ( x ) A(x) A(x)的次数的整数都是多项式 A ( x ) A(x) A(x)的一个次数界。
多项式加法
两个次数界为 n n n的多项式 A ( x ) A(x) A(x)和 B ( x ) B(x) B(x)的多项式相加,仍是一个次数界为 n n n的多项式 C ( x ) C(x) C(x),点值满足 C ( x ) = A ( x ) + B ( x ) C(x)=A(x)+B(x) C(x)=A(x)+B(x),且系数满足:
C ( x ) = ∑ j = 0 n − 1 ( a j + b j ) x j C(x) = \sum_{j=0}^{n-1} (a_j + b_j) x^j C(x)=j=0∑n−1(aj+bj)xj
多项式乘法
两个次数界为 n n n的多项式 A ( x ) A(x) A(x)和 B ( x ) B(x) B(x)的多项式相乘,其结果是一个次数界为 2 n − 1 2n-1 2n−1的多项式 C ( x ) C(x) C(x)。点值满足 C ( x ) = A ( x ) B ( x ) C(x)=A(x)B(x) C(x)=A(x)B(x),且系数满足:
C ( x ) = ∑ j = 0 2 n − 1 ∑ i = 0 j a i b j − i x j C(x) = \sum_{j=0}^{2n-1} \sum_{i=0}^ja_i b_{j-i} x^j C(x)=j=0∑2n−1i=0∑jaibj−ixj
也称卷积。
多项式表示法
系数表示
对于一个次数为 n n n的多项式 A ( x ) = ∑ j = 0 n − 1 a j x j A(x) = \sum_{j=0}^{n-1} a_j x^j A(x)=∑j=0n−1ajxj而言,其系数表示法是一组由系数组成的 n n n维列向量 a ⃗ = ( a 0 , a 1 , … , a n − 1 ) \vec{a} = (a_0,a_1,\ldots,a_{n-1}) a=(a0,a1,…,an−1)。
对于给定的一个系数表示法的多项式,我们可以通过霍纳法则(秦九韶算法)在 O ( n ) O(n) O(n)时间内计算出点值 A ( x 0 ) A(x_0) A(x0),即:
A ( x 0 ) = a 0 + x 0 ( a 1 + x 0 ( a 2 + … + x 0 ( a n − 1 ) ) … ) ) A(x_0) = a_0 + x_0 (a_1 +x_0(a_2+\ldots+x_0(a_{n-1}))\ldots)) A(x0)=a0+x0(a1+x0(a2+…+x0(an−1))…))
同时,我们对两个系数表示法的多项式只需要 O ( n ) O(n) O(n)的时间内就能完成,只需要做一次向量加法即可。
现在考虑系数表示法的多项式相乘,需要在 O ( n 2 ) O(n^2) O(n2)的时间内完成,因为需要做一次向量卷积(笛卡尔直积),记作 c ⃗ = a ⃗ ∗ b ⃗ \vec{c} = \vec{a} \ast \vec{b} c=a∗b。
点值表示
一个次数界为 n n n的多项式的点值表示就是由一个 n n n个点值对所组成的集合:
{ ( x 0 , y 0 ) , … , ( x n − 1 , y n − 1 ) } \{(x_0,y_0),\ldots,(x_{n-1},y_{n-1})\} {(x0,y0),…,(xn−1,yn−1)}
y k = A ( x k ) y_k = A(x_k) yk=A(xk),使得所有取值点的 x k x_k xk均不相同。
一个多项式的取值点的不同,表达也不同,故一个多项式的点值表示并不唯一。我们统一取值点为一个固定的点集 x 0 , … , x n − 1 x_0,\ldots,x_{n-1} x0,…,xn−1,称为基底。
求值计算的逆(从点值表示法计算系数表示法的过程)叫做插值。
定理:只有点值集合的大小和多项式的次数界相同,其插值才是唯一的。
证明:
点值和插值可以由如下矩阵表示:
[ 1 x 0 ⋯ x 0 n − 1 1 x 1 ⋯ x 1 n − 1 ⋮ ⋮ ⋱ ⋮ 1 x n − 1 ⋯ x n − 1 n − 1 ] [ a 0 a 1 ⋮ a n − 1 ] = [ y 0 y 1 ⋮ y n − 1 ] \begin{bmatrix} 1& x_{0}& \cdots & x_{0}^{n-1}\\ 1& x_{1}& \cdots & x_{1}^{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ 1& x_{n-1}& \cdots & x_{n-1}^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n-1 \end{bmatrix} =\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n-1 \end{bmatrix} ⎣ ⎡11⋮1x0x1⋮xn−1⋯⋯⋱⋯x0n−1x1n−1⋮xn−1n−1⎦ ⎤⎣ ⎡a0a1⋮an−1⎦ ⎤=⎣ ⎡y0y1⋮yn−1⎦ ⎤
左边的矩阵为范德蒙德矩阵,其行列式的值为 ∏ 0 ≤ j < k ≤ n − 1 ( x k − x j ) \prod_{0 \leq j \lt k \leq n-1}(x_k - x_j) ∏0≤j<k≤n−1(xk−xj)。
因为不存在两个相同的 x k x_k xk,因此该矩阵可逆,存在逆矩阵,该矩阵方程有唯一解。
该方法算作一个插值方法,通过LU等分解方式,加速矩阵乘积运算,可在 O ( n 3 ) O(n^3) O(n3)时间内进行插值。
作为了解,一种更快的插值方式是使用拉格朗日插值公式,其可以在 O ( n 2 ) O(n^2) O(n2)内完成插值。
因此,求值和插值是定义完备的互逆运算。
我们可以看到,通过点值计算具有相同基底多项式乘法的时间复杂度为 O ( n ) O(n) O(n),即基底不变 y k y_k yk对应相乘。同理,多项式加法的时间复杂度也为 O ( n ) O(n) O(n)。
计算一个点值的时间复杂度是 O ( n ) O(n) O(n)的,计算 n n n个不同的点的时间复杂度是 O ( n 2 ) O(n^2) O(n2)的,我们想法巧妙的构造一组基底,使得优化到 O ( n log n ) O(n \log n) O(nlogn)。
我们通过选取 2 n 2n 2n个单位复数根(记作 ω 2 n \omega_{2n} ω2n)作为基底,可以将求值和插值这两个过程优化到 O ( n log n ) O(n \log n) O(nlogn)。
其通过FFT计算多项式乘法定义了四个过程:
- 加倍次数界:因为结果的多项式的界是 2 n 2n 2n的,因此我们需要取 2 n 2n 2n个点。
- 求值:利用FFT算法,对两个多项式进行求值。
- 逐点相乘:对应点相乘。
- 插值:利用逆FFT算法,对多项式进行插值。
1和3步的时间为 O ( n ) O(n) O(n),24步的时间为 O ( n log n ) O(n \log n) O(nlogn)。
DFT与FFT
单位复数根
n次单位复数根是满足 ω n = 1 \omega^n=1 ωn=1的复数 ω \omega ω,这样的不同的复数恰好有 n n n个,对于 k = 0 , … , n − 1 k=0,\ldots,n-1 k=0,…,n−1, ω = e 2 k π n i \omega=e^{\frac{2k\pi}{n}i} ω=en2kπi。为了解释这个表达式,我们利用复数指数的形式解释:
e i u = cos ( u ) + i sin ( u ) e^{iu} = \cos(u) + i \sin(u) eiu=cos(u)+isin(u)
故 ω n = e 2 k π i = cos ( 2 k π ) + i sin ( 2 k π ) = 1 \omega^n=e^{2k\pi i}=\cos(2k\pi) + i\sin(2k\pi) = 1 ωn=e2kπi=cos(2kπ)+isin(2kπ)=1,其中指数 2 k π n \frac{2k\pi}{n} n2kπ,恰好把一个单位圆周分成 n n n分,正好对应着 n n n个不同的解,当 k = 1 k=1 k=1时, ω n = e 2 π i n \omega_n = e^{\frac{2\pi i}{n}} ωn=en2πi,称为主 n n n次单位复根,其他解都是主 n n n次单位复根的幂次,故下文表示 n n n个不同的根时,用 ω n 0 , ω n 1 , … , ω n n − 1 \omega_n^0,\omega_n^1,\ldots,\omega_n^{n-1} ωn0,ωn1,…,ωnn−1表示。
一般的,n次单位负数根在乘法意义上形成一个群,因为三角函数的周期性,n次单位复数根也存在周期性。
关于单位复数根的引理和证明
消去引理
对于任何的整数 n ≥ 0 , k ≥ 0 , d > 0 n \geq 0,k \geq 0,d\gt 0 n≥0,k≥0,d>0,有:
ω d n d k = e 2 d k π d n i = e 2 k π n i = ω n k \omega^{dk}_{dn} = e^{\frac{2 dk \pi}{dn}i} = e^{\frac{2 k \pi}{n}i} = \omega^{k}_{n} ωdndk=edn2dkπi=en2kπi=ωnk
推论:
ω n n 2 = ω 2 = − 1 \omega^{\frac{n}{2}}_n = \omega_2 = -1 ωn2n=ω2=−1
折半引理
如果 n > 0 n \gt 0 n>0是偶数,那么 n n n个 n n n次单位复数根的平方的集合,等同于 n / 2 n/2 n/2个 n / 2 n/2 n/2次单位复数根的的集合。
根据消去引理 ( ω n k ) 2 = ω n 2 k (\omega_n^k)^2=\omega^k_{\frac{n}{2}} (ωnk)2=ω2nk,如果对每个根都进行平方,那么每个不同的数正好出现两次,因为:
( ω n k + n / 2 ) 2 = ω n 2 k + n = ω n 2 k ω n n = ω n 2 k = ( ω n k ) 2 = ω n / 2 k (\omega_n^{k + n/2})^2 = \omega_n^{2k+n} = \omega_n^{2k}\omega_n^n = \omega_n^{2k} = (\omega_n^k)^2 = \omega_{n/2}^k (ωnk+n/2)2=ωn2k+n=ωn2kωnn=ωn2k=(ωnk)2=ωn/2k
因此, ( ω n k + n / 2 ) 2 (\omega_n^{k + n/2})^2 (ωnk+n/2)2和 ( ω n k ) 2 (\omega_n^k)^2 (ωnk)2和 ω n / 2 k \omega_{n/2}^k ωn/2k相同。
注意,折半引理对于分治FFT是十分有帮助的,因为他保证了子问题的父问题的规模的一半。
求和引理
对于任意的整数 n > 0 n \gt 0 n>0,和不能被 n n n整除的非负整数 k k k,存在:
∑ j = 0 n − 1 ( ω n k ) j = 0 \sum_{j=0}^{n-1}(\omega_{n}^{k})^j = 0 j=0∑n−1(ωnk)j=0
证明使用等比数列求和公式即可,这里不再证明。
DFT
我们希望在 ω n 0 , ω n 1 , … , ω n n − 1 \omega_n^0,\omega_n^1,\ldots,\omega_n^{n-1} ωn0,ωn1,…,ωnn−1处计算 A ( ω n k ) A(\omega_n^k) A(ωnk)的值,即在 n n n次单位复数根处进行求值。
如果多项式 A ( x ) = ∑ j = 0 n − 1 a j x j A(x) = \sum_{j=0}^{n-1}a_jx^j A(x)=∑j=0n−1ajxj,那么 y k = A ( ω n k ) = ∑ j = 0 n − 1 a j ω n k j y_k=A(\omega_n^k) = \sum_{j=0}^{n-1}a_j \omega_n^{kj} yk=A(ωnk)=∑j=0n−1ajωnkj。
那么向量 y ⃗ = ( y 0 , … , y n − 1 ) \vec{y}=(y_0,\ldots,y_{n-1}) y=(y0,…,yn−1)就是系数向量 a ⃗ = ( a 0 , … , a n − 1 ) \vec{a}=(a_0,\ldots,a_{n-1}) a=(a0,…,an−1)的离散傅里叶变换(DFT)。记作 y ⃗ = DFT n ( a ⃗ ) \vec{y} = \text{DFT}_n(\vec{a}) y=DFTn(a)。
FFT
这里介绍一种快速傅里叶变换(FFT)的算法,利用单位复数根的性质,我们就可在 O ( n log n ) O(n \log n) O(nlogn)的时间内计算出 DFT n ( a ⃗ ) \text{DFT}_n(\vec{a}) DFTn(a),通篇假设 n n n都是 2 2 2的整数次幂。
FFT利用了分治策略,分别采用偶数下标和奇数下标的系数构造两个新的多项式:
A 0 ( x ) = a 0 + a 2 x + a 4 x 2 + … + a n − 2 x n / 2 − 1 A 1 ( x ) = a 1 + a 3 x + a 5 x 2 + … + a n − 1 x n / 2 − 1 A_0(x) = a_0 + a_2x + a_4x^2 + \ldots + a_{n-2}x^{n/2 - 1} \\ A_1(x) = a_1 + a_3x + a_5x^2 + \ldots + a_{n-1}x^{n/2 - 1} \\ A0(x)=a0+a2x+a4x2+…+an−2xn/2−1A1(x)=a1+a3x+a5x2+…+an−1xn/2−1
其中 A 0 ( x ) A_0(x) A0(x)包含偶数下标, A 1 ( x ) A_1(x) A1(x)包含奇数下标,于是有:
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 A A在 ω n 0 , ω n 1 , … , ω n n − 1 \omega_n^0,\omega_n^1,\ldots,\omega_n^{n-1} ωn0,ωn1,…,ωnn−1处的值就转换为求界为:
- n / 2 n/2 n/2的两个多项式 A 0 A_0 A0和 A 1 A_1 A1在点 ( ω n 0 ) 2 , ( ω n 1 ) 2 , … , ( ω n n − 1 ) 2 (\omega_n^0)^2,(\omega_n^1)^2,\ldots,(\omega_n^{n-1})^2 (ωn0)2,(ωn1)2,…,(ωnn−1)2处的取值
- 根据 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 > 0 n \gt 0 n>0是偶数,那么 n n n个 n n n次单位复数根的平方的集合,等同于 n / 2 n/2 n/2个 n / 2 n/2 n/2次单位复数根的的集合。这样,问题就变成了两个界为 n / 2 n/2 n/2,单位复根数量为 n / 2 n/2 n/2的FFT问题,此时问题变成了两个子问题,但问题规模减半。
算法:
typedef complex<double> dcpx;
void recur_fft(dcpx A[], dcpx B[], int n, int inv)
{
if (n == 1)
return;
int hf = n >> 1;
for (int i = 0; i < hf; i++)
{
B[i] = A[i << 1];
B[i + hf] = A[(i << 1) | 1];
}
for (int i = 0; i < n; i++)
A[i] = B[i];
recur_fft(A, B, hf, inv);
recur_fft(A + hf, B + hf, hf, inv);
dcpx omega_n(cos(2 * M_PI / n), inv * sin(2 * M_PI / n));
dcpx omega(1, 0);
for (int i = 0; i < hf; i++)
{
B[i] = A[i] + omega * A[i + hf];
B[i + hf] = A[i] - omega * A[i + hf];
omega *= omega_n;
}
for (int i = 0; i < n; i++)
A[i] = B[i];
}
前几行按照奇偶性分类系数,然后递归,主要在最后循环合并答案上面,该代码 A A A代表输入的系数向量, B B B代表辅助数组,运行的结果将变换的结果放在 A A A向量里,参数 i n v inv inv稍后再说他的作用,在这里我们传入 1 1 1。
循环的最后一行,omega *= omega_n;
能保证我们能遍历
w
n
0
,
w
n
1
,
w
n
n
/
2
−
1
w_n^0,w_n^1,w_n^{n/2-1}
wn0,wn1,wnn/2−1。
考虑 k < n 2 k \lt \frac{n}{2} k<2n,此时表示成 k = i k = i k=i,那么 y 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 ) = y 0 [ i ] + ω n i y 1 [ i ] y_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_0[i] +\omega_n^i y_1[i] yk=A0(ωn2k)+ωnkA1(ωn2k)=A0(ωn/2k)+ωnkA1(ωn/2k)=y0[i]+ωniy1[i]。
考虑 k ≥ n 2 k \geq \frac{n}{2} k≥2n,此时表示成 k = i + n / 2 k = i + n/2 k=i+n/2,那么 y k = A 0 ( ω n / 2 i + n / 2 ) + ω n i + n / 2 A 1 ( ω n / 2 i + n / 2 ) = A 0 ( ω n / 2 i ) − ω n i A 1 ( ω n / 2 i ) = y 0 [ i ] − ω n i y 1 [ i ] y_k= A_0(\omega_{n/2}^{i+n/2}) + \omega_n^{i+n/2} A_1(\omega_{n/2}^{i+n/2}) = A_0(\omega_{n/2}^i) - \omega_n^i A_1(\omega_{n/2}^i) = y_0[i] - \omega_n^i y_1[i] yk=A0(ωn/2i+n/2)+ωni+n/2A1(ωn/2i+n/2)=A0(ωn/2i)−ωniA1(ωn/2i)=y0[i]−ωniy1[i]。
前面的系数 ± ω n i \pm \omega_n^i ±ωni,也称旋转因子。我们知道,乘以一个虚数相当于进行旋转。
根据主定理,其时间复杂度为 O ( n log n ) O(n \log n) O(nlogn)。
逆FFT
考虑复矩阵方程:
[ 1 ω 0 ⋯ ω 0 n − 1 1 ω 1 ⋯ ω 1 n − 1 ⋮ ⋮ ⋱ ⋮ 1 ω n − 1 ⋯ ω n − 1 n − 1 ] [ a 0 a 1 ⋮ a n − 1 ] = [ y 0 y 1 ⋮ y n − 1 ] \begin{bmatrix} 1& \omega_{0}& \cdots & \omega_{0}^{n-1}\\ 1& \omega_{1}& \cdots & \omega_{1}^{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ 1& \omega_{n-1}& \cdots & \omega_{n-1}^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n-1 \end{bmatrix} =\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n-1 \end{bmatrix} ⎣ ⎡11⋮1ω0ω1⋮ωn−1⋯⋯⋱⋯ω0n−1ω1n−1⋮ωn−1n−1⎦ ⎤⎣ ⎡a0a1⋮an−1⎦ ⎤=⎣ ⎡y0y1⋮yn−1⎦ ⎤
记最左边的矩阵为 V n V_n Vn,我们的目标就是找到 V n V_n Vn的逆矩阵 V n − 1 V_n^{-1} Vn−1,进而求解复矩阵方程。
定理:对于 j , k j,k j,k, V n [ j , k ] = ω n − j k n V_n[j,k] = \frac{\omega_n^{-jk}}{n} Vn[j,k]=nωn−jk。
证明:考虑证明 V n V n − 1 = I n V_nV_n^{-1}=I_n VnVn−1=In即可。需要使用上式求和引理。
给定逆矩阵 V n − 1 V_n^{-1} Vn−1,可以推导出 a j a_j aj:
a j = 1 n ∑ k = 0 n − 1 y k ω n − k j a_j = \frac{1}{n}\sum_{k=0}^{n-1}y_k \omega_n^{-kj} aj=n1k=0∑n−1ykωn−kj
我们和式子 y k = A ( ω n k ) = ∑ j = 0 n − 1 a j ω n k j y_k=A(\omega_n^k) = \sum_{j=0}^{n-1}a_j \omega_n^{kj} yk=A(ωnk)=∑j=0n−1ajωnkj进行对比,发现逆FFT过程只多了一个负号和除以 n n n,我们可以稍微修改一下FFT就可以完成逆FFT,即只需要把 i n v inv inv传入 − 1 -1 −1,并在最终除以 n n n即可。
离散傅里叶变换卷积定理:对于任意两个长度为 n n n的向量 a a a和 b b b,其中 n n n是 2 2 2的 n n n次幂,则
a ∗ b = DFT 2 n − 1 ( DFT 2 n ( a ) ⋅ DFT 2 n ( b ) ) a \ast b =\text{DFT}_{2n}^{-1}(\text{DFT}_{2n}(a) \cdot \text{DFT}_{2n}(b)) a∗b=DFT2n−1(DFT2n(a)⋅DFT2n(b))
其中向量 a a a和 b b b通过系数补零的方式扩大值 2 n 2n 2n即可。
代码实现
typedef complex<double> dcpx;
void recur_fft(dcpx A[], dcpx B[], int n, int inv)
{
if (n == 1)
return;
int hf = n >> 1;
for (int i = 0; i < hf; i++)
{
B[i] = A[i << 1];
B[i + hf] = A[(i << 1) | 1];
}
for (int i = 0; i < n; i++)
A[i] = B[i];
recur_fft(A, B, hf, inv);
recur_fft(A + hf, B + hf, hf, inv);
dcpx omega_n(cos(2 * M_PI / n), inv * sin(2 * M_PI / n));
dcpx omega(1, 0);
for (int i = 0; i < hf; i++)
{
B[i] = A[i] + omega * A[i + hf];
B[i + hf] = A[i] - omega * A[i + hf];
omega *= omega_n;
}
for (int i = 0; i < n; i++)
A[i] = B[i];
}
const int MAXT = 1 << 21;
dcpx A[MAXT];
dcpx B[MAXT];
dcpx T[MAXT];
int main()
{
FR;
int n, m;
scanf("%d %d", &n, &m);
int mx = max(n, m);
for (int i = 0; i < n + 1; i++)
{
int coff;
scanf("%d", &coff);
A[i].real(coff);
}
for (int i = 0; i < m + 1; i++)
{
int coff;
scanf("%d", &coff);
B[i].real(coff);
}
recur_fft(A, T, MAXT, 1);
recur_fft(B, T, MAXT, 1);
for (int i = 0; i < MAXT; i++)
{
A[i] *= B[i];
}
recur_fft(A, T, MAXT, -1);
for (int i = 0; i < n + m + 1; i++)
{
ll ans = A[i].real() / MAXT + 0.5;
printf("%lld ", ans);
}
return 0;
}
改进FFT
迭代版本的FFT
考虑我们的递归版的递归过程,我们是将其进行系数向量的拆解,然后逐层合并。如果我们可以在一开始就将系数向量按位置拆解,那么我们就可以通过循环的方法自底向上的方式进行合并答案。
考虑 n = 8 n=8 n=8的情况,重新排列系数下标依次为 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 a0,a4,a2,a6,a1,a5,a3,a7。我们看看 0 , 4 , 2 , 6 , 1 , 5 , 3 , 7 0,4,2,6,1,5,3,7 0,4,2,6,1,5,3,7和 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 0,1,2,3,4,5,6,7 0,1,2,3,4,5,6,7的关系,我们发现两个数列的二进制每一项都是对应项的转置。这种方法叫做逆位置对换。
我们可以通过递推的方式来实现:
for (int i = 1; i < MAXT; i++)
{
rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (E - 1));
}
考虑逆位置对换算法,我们遍历每一个位置 i i i,如果 i < r e v [ i ] i < rev[i] i<rev[i],说明 i i i应该交换,且一次交换正好能归位两个元素,因为 r e v [ r e v [ i ] ] = i rev[rev[i]] = i rev[rev[i]]=i。
for (int i = 0; i < n; i++)
if (i < rev[i])
swap(A[i], A[rev[i]]);
完整代码:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define FR freopen("in.txt", "r", stdin)
typedef long long ll;
typedef complex<double> dcpx;
const int E = 21;
const int MAXT = 1 << E;
const double PI = acos(-1);
int rev[MAXT];
void iter_fft(dcpx A[], int LOG, int inv)
{
int n = 1 << LOG;
for (int i = 0; i < n; i++)
if (i < rev[i])
swap(A[i], A[rev[i]]);
for (int e = 1; e <= LOG; e++)
{
int m = 1 << e;
for (int i = 0; i < n; i += m)
{
int hf = m / 2;
dcpx omega(1, 0);
dcpx omega_n(cos(2 * PI / m), inv * sin(2 * PI / m));
for (int j = 0; j < hf; j++)
{
dcpx x = A[i + j];
dcpx y = A[i + j + hf] * omega;
A[i + j] = x + y;
A[i + j + hf] = x - y;
omega *= omega_n;
}
}
}
}
dcpx A[MAXT];
dcpx B[MAXT];
int main()
{
for (int i = 1; i < MAXT; i++)
{
rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (E - 1));
}
int n, m;
cin >> n >> m;
for (int i = 0; i < n + 1; i++)
{
double a;
cin >> a;
A[i].real(a);
}
for (int i = 0; i < m + 1; i++)
{
double a;
cin >> a;
B[i].real(a);
}
iter_fft(A, E, 1);
iter_fft(B, E, 1);
for (int i = 0; i < MAXT; i++)
{
A[i] *= B[i];
}
iter_fft(A, E, -1);
for (int i = 0; i < m + n + 1; i++)
{
cout << ll(A[i].real() / MAXT + 0.5) << " ";
}
return 0;
}
例题
高精度乘法,考虑一个数字的按权展开的形式:
a = ∑ i = 0 … a i × 1 0 i a = \sum_{i=0}^{\dots}a_i \times 10^i a=i=0∑…ai×10i
这是一个多项式,那么两个数的乘积就是两个多项式相乘,此时就可以使用FFT卷起来了。
最后在考虑进位问题即可。
const int E = 21;
const int MAXT = 1 << E;
int rev[MAXT];
void iter_fft(dcpx A[], int inv)
{
for (int i = 0; i < MAXT; i++)
{
if (i < rev[i])
swap(A[i], A[rev[i]]);
}
for (int e = 1; e <= E; e++)
{
int m = 1 << e;
int hf = m >> 1;
for (int i = 0; i < MAXT; i += m)
{
dcpx omega(1, 0);
dcpx omega_n(cos(2 * M_PI / m), inv * sin(2 * M_PI / m));
for (int k = 0; k < hf; k++)
{
dcpx x = A[i + k];
dcpx y = A[i + k + hf] * omega;
A[i + k] = x + y;
A[i + k + hf] = x - y;
omega *= omega_n;
}
}
}
}
dcpx A[MAXT];
dcpx B[MAXT];
int main()
{
FR;
for (int i = 1; i < MAXT; i++)
{
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (E - 1));
}
string a, b;
cin >> a >> b;
for (int i = a.size() - 1; i >= 0; i--)
{
A[a.size() - i - 1] = a[i] - '0';
}
for (int i = b.size() - 1; i >= 0; i--)
{
B[b.size() - i - 1] = b[i] - '0';
}
iter_fft(A, 1);
iter_fft(B, 1);
for (int i = 0; i < MAXT; i++)
{
A[i] *= B[i];
}
iter_fft(A, -1);
string ans;
ll carry = 0;
for (int i = 0; i < MAXT; i++)
{
ll pp = ll(A[i].real() / MAXT + 0.5) + carry;
carry = pp / 10;
ans.push_back((pp % 10) + '0');
}
int i = MAXT - 1;
while (i >= 0 && ans[i] == '0')
i--;
string Ans;
for (; i >= 0; i--)
{
Ans.push_back(ans[i]);
}
cout << Ans;
return 0;
}
总结
本文介绍了多项式的基础和FFT算法的实现,是多项式算法系列的基础。