运算法则
算法
时间复杂度
多项式乘法
快速傅里叶变换
Θ ( n log 2 n ) \Theta(n\log_2 n) Θ(nlog2n)
倍增+快速数论变换
Θ ( n log 2 n ) \Theta(n\log_2 n) Θ(nlog2n)
求导+积分
Θ ( n log 2 n ) \Theta(n\log_2 n) Θ(nlog2n)
泰勒展开+牛顿迭代
Θ ( n log 2 n ) \Theta(n\log_2 n) Θ(nlog2n)
分治FFT/多项式求逆
Θ ( n log 2 2 n ) / Θ ( n log 2 n ) \Theta(n\log_2^2 n)/\Theta(n\log_2 n) Θ(nlog22n)/Θ(nlog2n)
% 能不能看得懂我不知道,但我相信至少高中水平是完全没问题的(我也才初中…… 高中了…… )
本章为多项式全家桶基础中的基础——多项式乘法。
缩写
全称
作用
时间复杂度
DFT
离散傅立叶变换
时频域转换
O ( n 2 ) O(n^2) O(n2)
FFT
快速傅立叶变换
时频域转换 ( ( (有精度误差 ) ) )
O ( 大常数 + n log 2 n ) O({\small\texttt{大常数}}+n\log_2n) O(大常数+nlog2n)
NTT/FNTT
模意义下的时频域转换
O ( 小常数 + n log 2 n ) O({\small\texttt{小常数}}+n\log_2n) O(小常数+nlog2n)
MTT
任意模数的NTT
任意模意义下的时频域转换
O ( n log 2 n ) O(n\log_2n) O(nlog2n)
FWT
快速沃尔什变换
快速集合卷积
O ( 不定 ) O({\small\texttt{不定}}) O(不定)
FMT
快速莫比乌斯变换
逆莫比乌斯反演?
O ( 不定 ) O({\small\texttt{不定}}) O(不定)
% 本文不包含但必不可少的前置技能:
% 什么是卷积?卷积是定义在函数上的运算。下面是百度百科的定义。
% 卷积是两个变量在某范围内相乘后求和的结果。如果卷积的变量是序列 x ( n ) x(n) x(n) 和 h ( n ) h(n) h(n),则卷积的结果为一个函数 y ( n ) = ∑ i = − ∞ ∞ x ( i ) h ( n − i ) = x ( n ) ∗ h ( n ) y(n)=\sum_{i=-\infty}^{\infty}x(i)h(n-i)=x(n)*h(n) y(n)=i=−∞∑∞x(i)h(n−i)=x(n)∗h(n)
% 出去的童鞋不要被百科拐跑了,百科哪有这里的文章没节操是不是…
修信号与系统的同学可能见过下面这段话:
% 给定义在 ( − ∞ , + ∞ ) (-\infty,+\infty) (−∞,+∞) 上的函数 f 1 ( t ) f_1(t) f1(t) 与 f 2 ( t ) f_2(t) f2(t) ,称由含参变量 t t t 的广义积分所确定的函数 g ( t ) = ∫ − ∞ + ∞ f 1 ( τ ) f 2 ( t − τ ) d τ g(t)=\int^{+\infty}_{-\infty}f_1(τ)f_2(t-τ)\ dτ g(t)=∫−∞+∞f1(τ)f2(t−τ) dτ 为函数 f 1 ( t ) f_1(t) f1(t) 与 f 2 ( t ) f_2(t) f2(t) 的卷积,记为: g ( t ) = f 1 ( t ) ∗ f 2 ( t ) g(t)=f_1(t)*f_2(t) g(t)=f1(t)∗f2(t)
% 这里最简单的一句话概括:多项式乘法实际上是多项式系数向量的卷积。还是很晕?举个栗子吧。
设有两个多项式 f = x 2 + 5 x + 4 , g = 3 x 2 − 5 x + 7 f=x^2+5x+4,g=3x^2-5x+7 f=x2+5x+4,g=3x2−5x+7。 则他们的卷积 h = f × g = 3 x 4 + 10 x 3 − 5 x 2 + 15 x + 28 h=f\times g=3x^4+10x^3-5x^2+15x+28 h=f×g=3x4+10x3−5x2+15x+28 其实是向量 f ⃗ = ( 1 , 5 , 4 ) , g ⃗ = ( 3 , − 5 , 7 ) \vec f=(1,5,4),\vec g=(3,-5,7) f =(1,5,4),g =(3,−5,7) 的卷积,为 h ⃗ = f ⃗ ⊗ g ⃗ = ( 3 , 10 , − 5 , 15 , 28 ) \vec h=\vec f\otimes \vec g=(3,10,-5,15,28) h =f ⊗g =(3,10,−5,15,28) 根据多项式乘法的计算方法,可得: h n = ∑ i = 0 n f i ⋅ g n − i h_n=\sum\limits_{i=0}^{n} f_i\cdot g_{n-i} hn=i=0∑nfi⋅gn−i 现在再回去看百度百科的词条(就看我节选的那一段),懂了吗?
% 设 f ( x ) f(x) f(x) 为一个 n n n 次函数,则其解析式为一个关于 x x x 的 n n n 次多项式。其解析式即为其系数表示法,也就是时域。
例: f ( x ) = x 2 + 5 x + 4 = ( 1 , 5 , 4 ) f(x)=x^2+5x+4=(1,5,4) f(x)=x2+5x+4=(1,5,4)
现在已知两个函数系数表示法,求这两个函数的卷积的系数表示法。显然地,需要 O ( n 2 ) O(n^2) O(n2) 的时间。
设有两个函数 f ( x ) = x 2 + 5 x + 4 , g ( x ) = 3 x 2 − 5 x + 7 f(x)=x^2+5x+4,g(x)=3x^2-5x+7 f(x)=x2+5x+4,g(x)=3x2−5x+7,则有
h ( x ) = f ( x ) ⊗ g ( x ) = ( x 2 + 5 x + 4 ) × ( 3 x 2 − 5 x + 7 ) = 3 x 4 + 10 x 3 − 5 x 2 + 15 x + 28 \begin{aligned}h(x)&=f(x)\otimes g(x)\\ &=(x^2+5x+4)\times (3x^2-5x+7)\\ &=3x^4+10x^3-5x^2+15x+28\end{aligned} h(x)=f(x)⊗g(x)=(x2+5x+4)×(3x2−5x+7)=3x4+10x3−5x2+15x+28
% 设 f ( x ) f(x) f(x) 为一个 n n n 次函数,因为 n + 1 n+1 n+1 个点确定一条 n n n 次函数。因此可以选取函数上的 n + 1 n+1 n+1 个点来表示一个函数,称为点值表示法,也就是频域。例如: f ( x ) = { ( − 1 , 0 ) , ( 0 , 4 ) , ( 1 , 10 ) } f(x)=\{(-1,0),(0,4),(1,10)\} f(x)={(−1,0),(0,4),(1,10)} 如果已知两个函数在相同 x x x 坐标上的点值表示法呢?求这两个函数的卷积的点值表示法。可以很负责任地说,只需要 O ( n ) O(n) O(n) 的时间就可以完成。令
f ( x ) = { ( − 1 , 0 ) , ( 0 , 4 ) , ( 1 , 10 ) } , g ( x ) = { ( − 1 , 15 ) , ( 0 , 7 ) , ( 1 , 5 ) } f(x)=\{(-1,0),(0,4),(1,10)\},g(x)=\{(-1,15),(0,7),(1,5)\} f(x)={(−1,0),(0,4),(1,10)},g(x)={(−1,15),(0,7),(1,5)} 则有
h ( x ) = f ( x ) ⊗ g ( x ) = { ( − 1 , 0 ) , ( 0 , 4 ) , ( 1 , 10 ) } ⊗ { ( − 1 , 15 ) , ( 0 , 7 ) , ( 1 , 5 ) } = { ( − 1 , 0 × 15 ) , ( 0 , 4 × 7 ) , ( 1 , 10 × 5 ) } = { ( − 1 , 0 ) , ( 0 , 28 ) , ( 1 , 50 ) } \begin{aligned}h(x)&=f(x)\otimes g(x)\\ &=\{(-1,0),(0,4),(1,10)\}\otimes \{(-1,15),(0,7),(1,5)\}\\ &=\{(-1,0\times15),(0,4\times 7),(1,10\times 5)\}\\ &=\{(-1,0),(0,28),(1,50)\}\end{aligned} h(x)=f(x)⊗g(x)={(−1,0),(0,4),(1,10)}⊗{(−1,15),(0,7),(1,5)}={(−1,0×15),(0,4×7),(1,10×5)}={(−1,0),(0,28),(1,50)} 可是不对啊,两个二次多项式的乘积应该是四次多项式啊,怎么能用三个点确定呢? 没错,不能用三个点确定,那我们就用更多个点来确定。
h ( x ) = f ( x ) ⊗ g ( x ) = { ( − 2 , 20 ) , ( − 1 , 0 ) , ( 0 , 4 ) , ( 1 , 10 ) , ( 2 , 18 ) } ⊗ { ( − 2 , 29 ) , ( − 1 , 15 ) , ( 0 , 7 ) , ( 1 , 5 ) , ( 2 , 9 ) } = { ( − 2 , 20 × 29 ) , ( − 1 , 0 × 15 ) , ( 0 , 4 × 7 ) , ( 1 , 10 × 5 ) , ( 2 , 18 × 9 ) } = { ( − 2 , 580 ) ( − 1 , 0 ) , ( 0 , 28 ) , ( 1 , 50 ) , ( 2 , 162 ) } \begin{aligned}h(x)&=f(x)\otimes g(x)\\ &=\{(-2,20),(-1,0),(0,4),(1,10),(2,18)\}\otimes \{(-2,29),(-1,15),(0,7),(1,5),(2,9)\}\\ &=\{(-2,20\times 29),(-1,0\times15),(0,4\times 7),(1,10\times 5),(2,18\times 9)\}\\ &=\{(-2,580)(-1,0),(0,28),(1,50),(2,162)\}\end{aligned} h(x)=f(x)⊗g(x)={(−2,20),(−1,0),(0,4),(1,10),(2,18)}⊗{(−2,29),(−1,15),(0,7),(1,5),(2,9)}={(−2,20×29),(−1,0×15),(0,4×7),(1,10×5),(2,18×9)}={(−2,580)(−1,0),(0,28),(1,50),(2,162)} 但是一般情况下很少使用点值表示法,因此在使用这种方法之前需要进行转换,那转换的代价是多少? T ( n ) = Θ ( n ) × Θ ( n ) = Θ ( n 2 ) T(n)=\Theta(n)\times \Theta(n)=\Theta(n^2) T(n)=Θ(n)×Θ(n)=Θ(n2) 似真似幻,漫天喜悦后,徒留一地空想。好在,我们都在阴沟里,但仍有人仰望星空。J.W.库利和T.W.图基于1965年发明了快速傅里叶变换 (fast Fourier transform)。
% 同时具有大小和方向的量。如下图中的向量OB记作 O B → \overrightarrow{OB} OB
% 定义虚数单位 i 2 = − 1 \text{i}^2=-1 i2=−1(不用斜体以示区别)。则对于 a , b ∈ R a,b\in \R a,b∈R(实数域) , b ≠ 0 ,b\not=0 ,b=0,形如 z = a + b i z=a+b\text{i} z=a+bi 的数称为虚数。其中 a a a 为实部, b b b 为虚部。特别地,当且仅当 a = 0 , b ≠ 0 a=0,b\not=0 a=0,b=0 时,我们称 z z z 为纯虚数。虚数和实数统称为复数。
在笛卡尔坐标系中,把 x x x 轴当做实轴,把 y y y 轴当做虚轴(单位长为 i \text{i} i),这样的坐标系叫做复平面。
这上图中的向量对应了复数 2 + 3 i 2+3\text{i} 2+3i 。
struct complex{
double x,y;//实部和虚部
complex (double xx=0,double yy=0) {
x=xx,y=yy;
}
};
复数的辐角 若复数 z z z 所表示的向量与 x x x 轴正半轴的夹角为 α \alpha α,则称 θ = α + 2 k π , k ∈ Z \theta=\alpha+2k\pi,k\in \Z θ=α+2kπ,k∈Z(整数域) 为复数 z z z 的辐角,特别地,当 θ ∈ ( − π , π ] \theta\in(-\pi,\pi] θ∈(−π,π] 时,我们称 θ \theta θ 为复数 z z z 的辐角主值,也称主辐角,记作 arg ( z ) \arg(z) arg(z)。如上图, 2 + 3 i 2+3i 2+3i 的辐角为 θ + 2 k π , k ∈ Z \theta+2k\pi,k\in \Z θ+2kπ,k∈Z,辐角主值为 θ \theta θ。
% 加法运算:实部和虚部分别相加。例:
( 3 + 5 i ) + ( 8 − 4 i ) = 11 + i (3+5\text{i})+(8-4\text{i})=11+\text{i} (3+5i)+(8−4i)=11+i 减法运算:实部和虚部分别相减。例:
( 3 + 5 i ) − ( 8 − 4 i ) = − 5 + 9 i (3+5\text{i})-(8-4\text{i})=-5+9\text{i} (3+5i)−(8−4i)=−5+9i 乘法运算:多项式乘法。例:
( 3 + 5 i ) ( 8 − 4 i ) = 3 × 8 − 3 × 4 i + 5 × 8 i − 5 i × 4 i = 24 − 12 i + 40 i − 20 i 2 = 44 + 28 i \begin{aligned}(3+5\text{i})(8-4\text{i})&=3\times8-3\times 4\text{i}+5\times 8\text{i}-5\text{i}\times 4\text{i}\\ &=24-12\text{i}+40\text{i}-20\text{i}^2\\ &=44+28\text{i}\end{aligned} (3+5i)(8−4i)=3×8−3×4i+5×8i−5i×4i=24−12i+40i−20i2=44+28i
complex operator+(complex a,complex b) {
return complex(a.x+b.x , a.y+b.y);
}
complex operator-(complex a,complex b) {
return complex(a.x-b.x , a.y-b.y);
}
complex operator*(complex a,complex b) {
return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);
}
复数的模(长) 对于虚数 z = a + b i z=a+b\text{i} z=a+bi,定义 z z z 的模(长)为 ∣ z ∣ = a 2 + b 2 |z|=\sqrt{a^2+b^2} ∣z∣=a2+b2 其几何意义是对应向量的模长,也是对应复平面上的点到原点的距离。
复数的三角形式 根据极坐标系的相关知识,在确定复数 z z z 的辐角(主值) θ \theta θ 与模长 r = ∣ z ∣ r=|z| r=∣z∣ 后复数有唯一解
z = r ( cos θ + i sin θ ) z=r(\cos\theta+\text{i}\sin\theta) z=r(cosθ+isinθ)这便是复数的三角表示,显然有
a = r cos θ , b = r s i n θ a=r\cos \theta,b=rsin\theta a=rcosθ,b=rsinθ 尽管复数的三角形式在加减法时显得十分笨拙,但复数的三角形式在乘除法方面上极具优势,一般地,对于两个复数 x = r 1 ( cos θ 1 + i sin θ 1 ) , y = r 2 ( cos θ 2 + i sin θ 2 ) x=r_1(\cos\theta_1+\text{i}\sin\theta_1),y=r_2(\cos\theta_2+\text{i}\sin\theta_2) x=r1(cosθ1+isinθ1),y=r2(cosθ2+isinθ2)根据三角恒等变换的知识,我们有 x × y = r 1 r 2 ( cos θ 1 + i sin θ 1 ) ( cos θ 2 + i sin θ 2 ) = r 1 r 2 [ ( cos θ 1 cos θ 2 − sin θ 1 sin θ 2 ) + i ( sin θ 1 cos θ 2 + cos θ 1 sin θ 2 ) ] = r 1 r 2 [ cos ( θ 1 + θ 2 ) + i sin ( θ 1 + θ 2 ) ] \begin{aligned} x\times y&=r_1r_2(\cos\theta_1+\text{i}\sin\theta_1)(\cos\theta_2+\text{i}\sin\theta_2)\\ &=r_1r_2[(\cos\theta_1\cos\theta_2-\sin\theta_1\sin\theta_2)+\text{i}(\sin\theta_1\cos\theta_2+\cos\theta_1\sin\theta_2)]\\ &=r_1r_2[\cos(\theta_1+\theta_2)+\text{i}\sin(\theta_1+\theta_2)] \end{aligned} x×y=r1r2(cosθ1+isinθ1)(cosθ2+isinθ2)=r1r2[(cosθ1cosθ2−sinθ1sinθ2)+i(sinθ1cosθ2+cosθ1sinθ2)]=r1r2[cos(θ1+θ2)+isin(θ1+θ2)]这说明两个复数相乘得到的复数模长等于两个复数模长之积,得到的复数辐角等于两个复数的辐角之和。特别地,复数 z n z^n zn 的辐角为 n arg ( z ) n\arg(z) narg(z),模长为 ∣ z ∣ n |z|^n ∣z∣n。
除法也同理,得到的复数主辐角相减,模长相除。
x y = r 1 r 2 ( cos θ 1 + i sin θ 1 ) ( cos θ 2 + i sin θ 2 ) = r 1 r 2 ( cos θ 1 + i sin θ 1 ) ( cos θ 2 − i sin θ 2 ) ( cos θ 2 + i sin θ 2 ) ( cos θ 2 − i sin θ 2 ) = r 1 r 2 [ ( cos θ 1 cos θ 2 + sin θ 1 sin θ 2 ) + i ( sin θ 1 cos θ 2 − cos θ 1 sin θ 2 ) cos 2 θ 2 + sin 2 θ 2 ] = r 1 r 2 [ cos ( θ 1 − θ 2 ) + i sin ( θ 1 − θ 2 ) ] \begin{aligned} \frac xy&=\frac {r_1}{r_2} \frac{(\cos\theta_1+\text{i}\sin\theta_1)}{(\cos\theta_2+\text{i}\sin\theta_2)}\\ &=\frac {r_1}{r_2} \frac{(\cos\theta_1+\text{i}\sin\theta_1)(\cos\theta_2-\text{i}\sin\theta_2)}{(\cos\theta_2+\text{i}\sin\theta_2)(\cos\theta_2-\text{i}\sin\theta_2)}\\ &=\frac {r_1}{r_2}[ \frac{(\cos\theta_1\cos\theta_2+\sin\theta_1\sin\theta_2)+\text{i}(\sin\theta_1\cos\theta_2-\cos\theta_1\sin\theta_2)}{\cos^2\theta_2+\sin^2\theta_2}]\\ &=\frac {r_1}{r_2}[\cos(\theta_1-\theta_2)+\text{i}\sin(\theta_1-\theta_2)] \end{aligned} yx=r2r1(cosθ2+isinθ2)(cosθ1+isinθ1)=r2r1(cosθ2+isinθ2)(cosθ2−isinθ2)(cosθ1+isinθ1)(cosθ2−isinθ2)=r2r1[cos2θ2+sin2θ2(cosθ1cosθ2+sinθ1sinθ2)+i(sinθ1cosθ2−cosθ1sinθ2)]=r2r1[cos(θ1−θ2)+isin(θ1−θ2)]
代数(学)基本定理 任何复系数一元 n n n 次多项式 方程在复数域上至少有一根( n ⩾ 1 n\geqslant 1 n⩾1)。对这个定理的证明目前无法绕开高等数学,由这个定理,我们可以发现任何实系数一元 n n n 次多项式方程在复数域上必有 n n n 个复数根,只需根据代数基本定理分解因式分解出所有因式即可。
% 一般地,关于 ω \omega ω 的方程 ω n = 1 , n ∈ Z \omega^n=1,n\in \Z ωn=1,n∈Z 在实数域内有且仅有一个根 1 1 1,由代数基本定理,在复数域内,这个方程有 n n n 个根,它们统称为 n n n 次单位根,本文用 ω n i \omega_n^i ωni 来表示最小非负辐角第 i + 1 i+1 i+1 大的 n n n 次单位根,可知 ω n 0 = 1 \omega_n^0=1 ωn0=1。
事实上,这个方程在复数域内的 n n n个根为 ω n i = cos ( i ⋅ 2 π n ) + i sin ( i ⋅ 2 π n ) \omega_n^i=\cos(i\cdot\frac {2\pi}{n})+\text{i}\sin(i\cdot\frac{2\pi}{n}) ωni=cos(i⋅n2π)+isin(i⋅n2π) ,因为其 n n n 次乘方对应的辐角为 n ⋅ i ⋅ 2 π n = 2 i π n\cdot i\cdot\frac {2\pi}{n}=2i\pi n⋅i⋅n2π=2iπ,故 arg ( ( ω n 1 ) n ) = 0 \arg\left((\omega_n^1)^n\right)=0 arg((ωn1)n)=0, ∣ ( ω n 1 ) n ∣ = 1 |(\omega_n^1)^n|=1 ∣(ωn1)n∣=1,因而 ( ω n 1 ) n = 1 (\omega_n^1)^n=1 (ωn1)n=1,满足方程。
把它们画在复平面上可以发现,将单位圆 n n n 等分,做 n n n 个向量,那么这 n n n 个向量对应复数即为 n n n 次单位根,以 1 1 1 为起点,逆时针分别为 ω n i \omega_n^i ωni。
例如当 n = 8 n=8 n=8 时,在复平面上表示为(圆为单位圆)
% 规定这样的次序还可以带来一些便利,例如 ω n i = ( ω n 1 ) i \omega_n^i=(\omega_n^1)^i ωni=(ωn1)i,不妨再规定 ω n 1 = ω n \omega_n^1=\omega_n ωn1=ωn,那么便可以使得 ω n i \omega_n^i ωni 既是第 i i i 个 n n n 次单位根,又表示 ω n \omega_n ωn 的 i i i 次方,同时将 ω n \omega_n ωn 的上标由 [ 0 , n − 1 ] [0,n-1] [0,n−1] 中的整数扩充到了 Z \Z Z。
性质1 w n 0 = w n n = 1 w^0_n=w_n^n=1 wn0=wnn=1
性质2 w n 0 , w n 1 , ⋯ , w n n − 1 w_n^0, w_n^1, \cdots,w_n^{n-1} wn0,wn1,⋯,wnn−1 互不相同。
性质3 w 2 n 2 k = w n k w^{2k}_{2n}=w^k_n w2n2k=wnk
证 代入定义式即可,也可直接观察去掉单位圆上的奇数次单位根。
w 2 n 2 k = cos 2 k × 2 π 2 n + i sin 2 k × 2 π 2 n = cos k × 2 π n + i sin k × 2 π n = w n k \begin{aligned} w_{2n}^{2k}=&\cos 2k\times \frac{2\pi}{2n}+\text{i}\sin 2k\times \frac{2\pi}{2n}\\ =&\cos k\times \dfrac{2\pi}n+\text{i}\sin k\times \dfrac{2\pi}n\\ =&\ w_n^k\\ \end{aligned} w2n2k===cos2k×2n2π+isin2k×2n2πcosk×n2π+isink×n2π wnk性质4 w n k + n 2 = − w n k w_n^{k+\frac{n}{2}}=-w_n^k wnk+2n=−wnk
证 代入定义式即可,也可看出旋转半圈前后关于原点对称。
w n k + n 2 = w n k × w n n 2 = w n k × [ cos ( n 2 × 2 π n ) + i sin ( n 2 × 2 π n ) ] = w n k × ( cos π + i sin π ) = w n k × ( − 1 + 0 ) = − w n k \begin{aligned} w_n^{k+\frac{n}{2}}=&\ w_n^k\times w_n^{\frac{n}{2}}\\ =&\ w_n^k\times \left[\cos\left(\frac n2\times \frac{2\pi}{n}\right)+\text{i}\sin\left(\dfrac n2\times\dfrac{2\pi}{n}\right)\right]\\ =&\ w_n^k\times (\cos\pi+\text{i}\sin\pi)\\ =&\ w_n^k\times (-1+0)\\ =&-w_n^k \end{aligned} wnk+2n===== wnk×wn2n wnk×[cos(2n×n2π)+isin(2n×n2π)] wnk×(cosπ+isinπ) wnk×(−1+0)−wnk性质5 ∑ i = 0 n − 1 ( ω n j − k ) i = { 0 , k ≠ j n , k = j \begin{aligned}\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i =\begin{cases}0,k\ne j\\n,k=j\\\end{cases}\end{aligned} i=0∑n−1(ωnj−k)i={0,k=jn,k=j
证 1. 当 k ≠ j k\ne j k=j 时,根据等比数列的求和公式,可得:
∑ i = 0 n − 1 ( ω n j − k ) i = ( w n j − k ) n − 1 w n j − k − 1 = ( w n n ) j − k − 1 w n j − k − 1 = 1 − 1 w n j − k − 1 = 0 \begin{aligned} \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i =&\dfrac{(w_n^{j-k})^{n}-1}{w_n^{j-k}-1} =\dfrac{(w_n^n)^{j-k}-1}{w_n^{j-k}-1} =\dfrac{1-1}{w_n^{j-k}-1} =\ 0\\\end{aligned} i=0∑n−1(ωnj−k)i=wnj−k−1(wnj−k)n−1=wnj−k−1(wnn)j−k−1=wnj−k−11−1= 0 2.当 k = j k=j k=j 时,可得: ∑ i = 0 n − 1 ( ω n j − k ) i = ∑ i = 0 n − 1 1 = n \begin{aligned}\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i =\sum\limits_{i=0}^{n-1}\ 1=n\\ \end{aligned} i=0∑n−1(ωnj−k)i=i=0∑n−1 1=n
某些奇奇怪怪的读法。
% FFT 最常见的算法是 Cooley-Tukey \texttt{Cooley-Tukey} Cooley-Tukey 算法,它的基本思路在 1965 年由 J.W.Cooley \texttt{J.W.Cooley} J.W.Cooley 和 J.W.Tukey \texttt{J.W.Tukey} J.W.Tukey 提出的,它是一个基于分治策略的算法。
% 显然地,一个 n n n 次多项式可以被 n + 1 n+1 n+1 个点唯一确定。很恶心地, Cooley-Tukey \texttt{Cooley-Tukey} Cooley-Tukey 算法在取点的时候不取整数点,不取有理数点,甚至还不取实数点!而是取复数点, n n n 次单位根的 0 0 0 到 n − 1 n-1 n−1 次幂。即使是这样,那时间复杂度仍然是 Θ ( n 2 ) \Theta(n^2) Θ(n2) 的啊,而且复数的运算还比整数慢!别着急,来推推柿子。
设函数 A ( x ) A(x) A(x) 的系数为 ( a 0 , a 1 , a 2 , . . . , a n − 1 ) (a_0,a_1,a_2,…,a_{n-1}) (a0,a1,a2,…,an−1)。则有: A ( x ) = ∑ i = 0 n − 1 a i x i A(x)=\sum\limits_{i=0}^{n-1} a_ix^i A(x)=i=0∑n−1aixi 按照下标奇偶性分类,则有
A ( x ) = ( a 0 + a 2 x 2 + a 4 x 4 + . . . + a n − 2 x n − 2 ) + ( a 1 x + a 3 x 3 + a 5 x 5 + . . . + a n − 1 x n − 1 ) A(x)=(a_0+a_2x^2+a_4x^4+…+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+…+a_{n-1}x^{n-1}) A(x)=(a0+a2x2+a4x4+…+an−2xn−2)+(a1x+a3x3+a5x5+…+an−1xn−1)
设 A 1 ( x ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n 2 − 1 A 2 ( x ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n 2 − 1 \begin{aligned}A_1(x)=a_0+a_2x+a_4x^2+…+a_{n-2}x^{\frac n2-1}\\ A_2(x)=a_1+a_3x+a_5x^2+…+a_{n-1}x^{\frac n2-1}\end{aligned} A1(x)=a0+a2x+a4x2+…+an−2x2n−1A2(x)=a1+a3x+a5x2+…+an−1x2n−1 则有: A ( x ) = A 1 ( x 2 ) + x A 2 ( x 2 ) A(x)=A_1(x^2)+xA_2(x^2) A(x)=A1(x2)+xA2(x2)
% 这个时候,单位根终于派上用场了!
没错,我们把 x = w n k ( k < n 2 ) x=w_n^k (k<\frac{n}{2}) x=wnk(k<2n) 代入,得
A ( w n k ) = A 1 ( w n 2 k ) + w n k A 2 ( w n 2 k ) = A 1 ( w n 2 k ) + w n k A 2 ( w n 2 k ) \begin{aligned}A(w_n^k)&=A_1(w_n^{2k})+w_n^kA_2(w_n^{2k})\\ &=A_1(w_{\frac n2}^{k})+w_{n}^kA_2(w_{\frac n2}^{k})\end{aligned} A(wnk)=A1(wn2k)+wnkA2(wn2k)=A1(w2nk)+wnkA2(w2nk)
% 将 x = w n k + n 2 x=w_n^{k+\frac{n}{2}} x=wnk+2n 代入,得:
A ( w n k + n 2 ) = A 1 ( w n 2 k + n ) + w n k + n 2 A 2 ( w n 2 k + n ) = A 1 ( w n 2 k ∗ w n n ) − w n k A 2 ( w n 2 k ∗ w n n ) = A 1 ( w n 2 k ) − w n k A 2 ( w n 2 k ) = A 1 ( w n 2 k ) − w n k A 2 ( w n 2 k ) \begin{aligned}A(w_n^{k+\frac{n}{2}})&=A_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}A_2(w_n^{2k+n})\\ &=A_1(w_n^{2k}*w_n^n)-w_n^kA_2(w_n^{2k}*w_n^n)\\ &=A_1(w_n^{2k})-w_n^kA_2(w_n^{2k})\\ &=A_1(w_{\frac n2}^{k})-w_{n}^kA_2(w_{\frac n2}^{k})\end{aligned} A(wnk+2n)=A1(wn2k+n)+wnk+2nA2(wn2k+n)=A1(wn2k∗wnn)−wnkA2(wn2k∗wnn)=A1(wn2k)−wnkA2(wn2k)=A1(w2nk)−wnkA2(w2nk) 发现了什么?只相差一个符号!有什么作用?那么当我们在枚举第一个式子的时候,我们可以 O ( 1 ) O(1) O(1) 得到第二个式子的值。而且第一个式子的 k k k 在取遍了 [ 0 , n 2 − 1 ] [0,\frac n2-1] [0,2n−1] 的时候,第二个式子取遍了 [ n 2 , n − 1 ] [\frac n2,n-1] [2n,n−1]。因此原问题的规模缩小了一半,然后呢?
分治就好啦!下面是代码:
void fast_fast_tle(int limit,complex *a) {
if(limit==1) //只有一个常数项
return ;
complex a1[limit>>1],a2[limit>>1];
for(int i=0; i<=limit; i+=2) //根据下标的奇偶性分类
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fast_fast_tle(limit>>1,a1);//分治
fast_fast_tle(limit>>1,a2);
complex Wn=complex(std::cos(2.0*Pi/limit),std::sin(2.0*Pi/limit)),w=complex(1,0);
//Wn为单位根,w表示幂
for(int i=0; i<(limit>>1); i++,w=w*Wn) //这里的w相当于公式中的k
a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];//利用单位根的性质,O(1)得到另一部分
}
% 时间复杂度: Θ ( 大常数 + n l o g 2 n ) \Theta({\small \texttt{大常数}}+nlog_2n) Θ(大常数+nlog2n)这个大常数是哪里来的呢?复数运算!复数的乘法是非常慢的。因此如果不是数据非常大,应尽量避免使用FFT。
% 经过了 Cooley-Tukey \texttt{Cooley-Tukey} Cooley-Tukey 算法的折磨之后呢?然后直接乘就好了啊。
fast_fast_tle(n,a);
fast_fast_tle(n,b);
for(int i=0;i<=n;i++)
a[i]=a[i]*b[i];
% 时间复杂度 Θ ( n ) \Theta(n) Θ(n)。完了?当然没有。还要逆回去。
% 设 ( y 0 , y 1 , . . . , y n − 1 ) (y_0,y_1,…,y_{n-1}) (y0,y1,…,yn−1) 是 ( a 0 , a 1 , a 2 , … , a n − 1 ) (a_0,a_1,a_2,\dots,a_{n-1}) (a0,a1,a2,…,an−1) 的傅里叶变换,即 ( y 0 , y 1 , . . . , y n − 1 ) (y_0,y_1,…,y_{n-1}) (y0,y1,…,yn−1) 是 ( a 0 , a 1 , a 2 , … , a n − 1 ) (a_0,a_1,a_2,\dots,a_{n-1}) (a0,a1,a2,…,an−1) 在 ( w n 0 , w n 1 , … , w n n − 1 ) (w_n^0,w_n^1,\dots,w_n^{n-1}) (wn0,wn1,…,wnn−1)处的值。则有
y i = a 0 + a 1 w n 1 + a 2 ( w n 1 ) 2 + ⋯ + a n ( w n 1 ) n = a 0 + a 1 w n 1 + a 2 w n 2 + ⋯ + a n w n n = ∑ j = 0 n − 1 a j ( w n i ) j \begin{aligned}y_i&=a_0+a_1w_n^1+a_2(w_n^1)^2+\dots +a_n(w_n^1)^n\\ &=a_0+a_1w_n^1+a_2w_n^2+\dots +a_nw_n^n\\ &=\sum\limits_{j=0}^{n-1}a_j(w_n^i)^j\end{aligned} yi=a0+a1wn1+a2(wn1)2+⋯+an(wn1)n=a0+a1wn1+a2wn2+⋯+anwnn=j=0∑n−1aj(wni)j 设 ( c 0 , c 1 , c 2 , … , c n − 1 ) (c_0,c_1,c_2,\dots,c_{n-1}) (c0,c1,c2,…,cn−1) 是 ( y 0 , y 1 , . . . , y n − 1 ) (y_0,y_1,…,y_{n-1}) (y0,y1,…,yn−1) 在 ( w n 0 , w n − 1 , … , w n − ( n − 1 ) ) (w_n^0,w_n^{-1},\dots,w_n^{-(n-1)}) (wn0,wn−1,…,wn−(n−1)) 处的取值。
% 求C序列同样用 Cooley-Tukey \texttt{Cooley-Tukey} Cooley-Tukey 算法实现,只需要将只需要单位根变为 w n − 1 = w n n − 1 w_n^{-1}=w_n^{n-1} wn−1=wnn−1,相当于顺时针旋转。可以发现,在复平面上, w n n − 1 w_n^{n-1} wnn−1 和 w n 1 w_n^1 wn1 的 x x x 坐标相同, y y y 坐标互为相反数,因此在代码中只需要把上面第9行的sin(2.0*Pi/limit)
改成-sin(2.0*Pi/limit)
就可以了。
c k = ∑ i = 0 n − 1 y i ( w n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( w n i ) j ) ( w n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( w n j ) i ) ( w n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( w n j ) i ( w n − k ) i ) = ∑ i = 0 n − 1 ∑ j = 0 n − 1 a j ( w n j ) i ( w n − k ) i = ∑ i = 0 n − 1 ∑ j = 0 n − 1 a j ( w n j − k ) i = ∑ j = 0 n − 1 a j ( ∑ i = 0 n − 1 ( w n j − k ) i ) ∴ c k = ∑ j = 0 n − 1 a j ( ∑ i = 0 n − 1 ( w n j − k ) i ) \begin{aligned} c_k&=\sum_{i=0}^{n-1}y_i(w_n^{-k})^i\\ &=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(w_n^i)^j)(w_n^{-k})^i\\ &=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(w_n^j)^i)(w_n^{-k})^i\\ &=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(w_n^j)^i(w_n^{-k})^i)\\ &=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j(w_n^j)^i(w_n^{-k})^i\\ &=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j(w_n^{j-k})^i\\ &=\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(w_n^{j-k})^i)\\ \therefore c_k&=\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(w_n^{j-k})^i)\end{aligned} ck∴ck=i=0∑n−1yi(wn−k)i=i=0∑n−1(j=0∑n−1aj(wni)j)(wn−k)i=i=0∑n−1(j=0∑n−1aj(wnj)i)(wn−k)i=i=0∑n−1(j=0∑n−1aj(wnj)i(wn−k)i)=i=0∑n−1j=0∑n−1aj(wnj)i(wn−k)i=i=0∑n−1j=0∑n−1aj(wnj−k)i=j=0∑n−1aj(i=0∑n−1(wnj−k)i)=j=0∑n−1aj(i=0∑n−1(wnj−k)i) 根据性质5,有:
% 因而有
c k n = a k \dfrac{c_k}n=a_k nck=ak 然后呢?做完啦!先用 特别的 Cooley-Tukey \texttt{Cooley-Tukey} Cooley-Tukey 算法求出 ( c 0 , c 1 , c 2 , … , c n − 1 ) (c_0,c_1,c_2,\dots,c_{n-1}) (c0,c1,c2,…,cn−1),然后再除以 n n n 就好啦!两种FFT代码有很多类似的部分,因此可以弄多一个参数,表示做的是那一种FFT。
void fast_fast_tle(int limit,complex *a,int type){
if(limit==1)
return ;
complex a1[limit>>1],a2[limit>>1];
for(int i=0;i<=limit;i+=2)
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fast_fast_tle(limit>>1,a1,type);
fast_fast_tle(limit>>1,a2,type);
//当type=1时,表示逆时针的FFT;当type=-1时,表示顺时针的FFT。
complex Wn=complex(std::cos(2.0*Pi/limit),type*std::sin(2.0*Pi/limit)),w=complex(1,0);
for(int i=0;i<(limit>>1);i++,w=w*Wn)
a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];
}
% 执行完后还要除以 n n n,即:
for(int i=1;i<=limit;i++)
a[i].x/=limit;
% 令 N N N 为大于 n + m n+m n+m 的最小的 2 2 2 的正整数次幂。
不难看出,FFT的时空复杂度均为 O ( 不小的常数 + N l o g 2 N ) O({\small\texttt{不小的常数}}+Nlog_2N) O(不小的常数+Nlog2N)
下面是完整的Luogu P3803代码
#include<cmath>
#include<cstdio>
const int MAXN=2*1e6+10;
const double Pi=std::acos(-1.0);
struct complex {
double x,y;
complex (double xx=0,double yy=0) {
x=xx,y=yy;
}
} a[MAXN],b[MAXN];
complex operator+(complex a,complex b) {
return complex(a.x+b.x,a.y+b.y);
}
complex operator-(complex a,complex b) {
return complex(a.x-b.x,a.y-b.y);
}
complex operator*(complex a,complex b) {
return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
void fast_fast_tle(int limit,complex *a,int type){
if(limit==1)
return ;
complex a1[limit>>1],a2[limit>>1];
for(int i=0;i<=limit;i+=2)
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fast_fast_tle(limit>>1,a1,type);
fast_fast_tle(limit>>1,a2,type);
complex Wn=complex(std::cos(2.0*Pi/limit),type*std::sin(2.0*Pi/limit)),w=complex(1,0);
for(int i=0;i<(limit>>1);i++,w=w*Wn)
a[i]=a1[i]+w*a2[i],a[i+(limit>>1)]=a1[i]-w*a2[i];
}
int main() {
int n,m;
scanf("%d%d",&n,&m);
for(int i=0; i<=n; i++)
scanf("%lf",&a[i].x);
for(int i=0; i<=m; i++)
scanf("%lf",&b[i].x);
int limit=1;
while(limit<=n+m) limit<<=1;
fast_fast_tle(limit,a,1);
fast_fast_tle(limit,b,1);
for(int i=0; i<=limit; i++)
a[i]=a[i]*b[i];
fast_fast_tle(limit,a,-1);
for(int i=0; i<=n+m; i++)
printf("%d ",(int)floor(a[i].x/limit+0.5));
//由于FFT涉及到浮点数运算,因此需要考虑精度误差。
return 0;
}
% 于是乎,Luogu的提交结果无外乎有两种,一种是RE,另一种也是RE。
注意到空间复杂度中有 O ( N l o g 2 N ) O(Nlog_2N) O(Nlog2N)的空间是开在栈空间里的,因此可能会爆栈!所以有了下一小节。
% 观察一下原下标和新下标的二进制,发现了什么规律?
规律不是很好阐述,这里定义 x x x 的 l l l 位翻转为 x x x 在长度为 log 2 l − 1 \log_2l-1 log2l−1 的二进制下的翻转1。
举个例子: 1 1 1 的 8 8 8 位翻转为 4 4 4。
1 1 1 的长度为 log 2 8 − 1 = 3 \log_28-1=3 log28−1=3 的二进制为 001 001 001。前后翻转为 100 100 100,十进制下即为 4 4 4。
可以发现,原下标为 x x x 的新下标为 x x x 的 n n n 位翻转。设其为 r [ x ] r[x] r[x]。
若已经求出了 r [ i ] r[i] r[i],则可以写出如下代码:
void fast_fast_tle(complex *A,int type) {
for(int i=0;i<limit;i++) //求出要迭代的序列
if(i<r[i]) swap(A[i],A[r[i]]);
for(int mid=1;mid<limit;mid<<=1) { //待合并区间的长度(块的长度)
complex Wn(std::cos(Pi/mid),type*std::sin(Pi/mid));//单位根
for(int R=mid<<1,j=0;j<limit;j+=R) { //R是区间的右端点,j表示前已经到哪个位置了(枚举那一块)
complex w(1,0);//幂
for(int k=0;k<mid;k++,w=w*Wn) { //块内枚举
complex x=A[j+k],y=w*A[j+mid+k];
A[j+k]=x+y;
A[j+mid+k]=x-y;
}
}
}
}
% 那现在只剩下一个问题了: r [ i ] r[i] r[i] 怎么求?
% 由于FFT的时间复杂度为 O ( N l o g 2 N ) O(Nlog_2N) O(Nlog2N)。
% 因此如果求翻转序列的复杂度比 O ( N l o g 2 N ) O(Nlog_2N) O(Nlog2N) 还要高那就完了。
% 时间复杂度: O ( N l o g 2 N ) O(Nlog_2N) O(Nlog2N)
% 空间复杂度: O ( n ) O(n) O(n)
% 代码复杂度:小
% 代码:下面的 l = log 2 l i m i t l=\log_2 limit l=log2limit,实现上可以在求 l i m i t limit limit 的同时求出。
for(int i=1;i<=n-1;i++){
int t=i,ret=0;
for(int j=1;j<=l;j++){
ret=(ret<<1)|(t&1);
t>>=1;
}
r[i]=ret;
}
% 想不到吧,这种东西还有线性算法。
% 下面以 184 184 184 的 256 256 256 位翻转为例子。
% 我们把一个数的二进制分成两部分,设这个数的二进制长度为 l l l 位。
% 第一部分是前 l − 1 l-1 l−1 位,第二部分是最后一位。
1 0 1 1 1 0 0 ∣ 0 \qquad1\quad0\quad1\quad1\quad1\quad0\quad0\quad|\quad0 1011100∣0
% 184 184 184 的 256 256 256 位翻转等价于第二部分接上第一部分的 128 128 128 位翻转。即
0 ∣ 0 0 1 1 1 0 1 \qquad0\quad|\quad0\quad0\quad1\quad1\quad1\quad0\quad1 0∣0011101
% 那第一部分的 128 128 128 位翻转怎么求呢?等价于第一部分的 256 256 256 位翻转右移 1 1 1 位。
% 第一部分的 256 256 256 位翻转
0 0 1 1 1 0 1 0 \qquad0\quad0\quad1\quad1\quad1\quad0\quad1\quad0 00111010
% 右移一位,得:
0 0 1 1 1 0 1 \;\;\quad\qquad0\quad0\quad1\quad1\quad1\quad0\quad1 0011101
% 代码:
for(int i=1;i<limit;i++)
r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
时间复杂度: O ( n ) O(n) O(n)
空间复杂度: O ( n ) O(n) O(n)
代码复杂度:更小
Accepted -O2 \text{Accepted -O2} Accepted -O2 / 用时: 5111 m s 5111ms 5111ms / 内存: 243752 KB 243752\text{KB} 243752KB
#include<bits/stdc++.h>
using namespace std;
const int maxn=2148576;
const double pi=std::acos(-1.0);
struct comp{
double x,y;
comp(double xx=0,double yy=0):x(xx),y(yy) {}
friend comp operator+(const comp &x,const comp &y) {return comp(x.x+y.x,x.y+y.y);}
friend comp operator-(const comp &x,const comp &y) {return comp(x.x-y.x,x.y-y.y);}
friend comp operator*(const comp &a,const comp &b) {return comp(a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y);}
}a[maxn],b[maxn];
int limit=1,n,m,l=0,r[maxn];
void fft(comp *t,int ty){
for(int i=0;i<limit;i++)
if(i<r[i])
swap(t[i],t[r[i]]);
for(int mid=1;mid<limit;mid<<=1){
comp wn(std::cos(pi/mid),ty*std::sin(pi/mid));
for(int j=0,R=(mid<<1);j<limit;j+=R){
comp w(1,0);
for(int k=0;k<mid;k++,w=w*wn){
comp x=t[j+k],y=w*t[j+k+mid];
t[j+k]=x+y;
t[j+k+mid]=x-y;
}
}
}
}
int main(void)
{
scanf("%d%d",&n,&m);
while(limit<=n+m)
limit<<=1,++l;
for(int i=1;i<limit;i++)
r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
for(int i=0;i<=n;++i)
scanf("%lf",&a[i].x);
for(int i=0;i<=m;++i)
scanf("%lf",&b[i].x);
fft(a,1);
fft(b,1);
for(int i=0;i<limit;i++)
a[i]=a[i]*b[i];
fft(a,-1);
for(int i=0;i<=n+m;i++)
printf("%d ",(int)(a[i].x/limit+0.5));
return 0;
}
% 恭喜完成一万字阅读。
题目
来源
题解
A × \times ×B Problem
[MUTC2013]idiots
手机扫一扫
移动阅读更方便
你可能感兴趣的文章