快速傅里叶变换
阅读原文时间:2022年02月08日阅读:1

多项式全家桶

运算法则

算法

时间复杂度

多项式乘法

快速傅里叶变换

Θ ( n log ⁡ 2 n ) \Theta(n\log_2 n) Θ(nlog2​n)

多项式求逆

倍增+快速数论变换

Θ ( n log ⁡ 2 n ) \Theta(n\log_2 n) Θ(nlog2​n)

多项式对数函数

求导+积分

Θ ( n log ⁡ 2 n ) \Theta(n\log_2 n) Θ(nlog2​n)

多项式指数函数

泰勒展开+牛顿迭代

Θ ( n log ⁡ 2 n ) \Theta(n\log_2 n) Θ(nlog2​n)

分治FFT卷积

分治FFT/多项式求逆

Θ ( n log ⁡ 2 2 n ) / Θ ( n log ⁡ 2 n ) \Theta(n\log_2^2 n)/\Theta(n\log_2 n) Θ(nlog22​n)/Θ(nlog2​n)


作者的话

%   能不能看得懂我不知道,但我相信至少高中水平是完全没问题的(我也才初中…… 高中了…… )
  本章为多项式全家桶基础中的基础——多项式乘法。


分类

缩写

全称

作用

时间复杂度

DFT

离散傅立叶变换

时频域转换

O ( n 2 ) O(n^2) O(n2)

FFT

快速傅立叶变换

时频域转换 ( ( (有精度误差 ) ) )

O ( 大常数 + n log ⁡ 2 n ) O({\small\texttt{大常数}}+n\log_2n) O(大常数+nlog2​n)

NTT/FNTT

快速数论变换

模意义下的时频域转换

O ( 小常数 + n log ⁡ 2 n ) O({\small\texttt{小常数}}+n\log_2n) O(小常数+nlog2​n)

MTT

任意模数的NTT

任意模意义下的时频域转换

O ( n log ⁡ 2 n ) O(n\log_2n) O(nlog2​n)

FWT

快速沃尔什变换

快速集合卷积

O ( 不定 ) O({\small\texttt{不定}}) O(不定)

FMT

快速莫比乌斯变换

逆莫比乌斯反演?

O ( 不定 ) O({\small\texttt{不定}}) O(不定)


快速傅里叶变换

前置技能

%   本文不包含但必不可少的前置技能:

  1. sin函数,cos函数,弧度制
  2. 函数
  3. 笛卡尔坐标系(也叫做平面直角坐标系)
  4. ∑ \sum ∑ 的性质

卷积

%   什么是卷积?卷积是定义在函数上的运算。下面是百度百科的定义。

%   卷积是两个变量在某范围内相乘后求和的结果。如果卷积的变量是序列 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∑n​fi​⋅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​=r1​r2​(cosθ1​+isinθ1​)(cosθ2​+isinθ2​)=r1​r2​[(cosθ1​cosθ2​−sinθ1​sinθ2​)+i(sinθ1​cosθ2​+cosθ1​sinθ2​)]=r1​r2​[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​​=r2​r1​​(cosθ2​+isinθ2​)(cosθ1​+isinθ1​)​=r2​r1​​(cosθ2​+isinθ2​)(cosθ2​−isinθ2​)(cosθ1​+isinθ1​)(cosθ2​−isinθ2​)​=r2​r1​​[cos2θ2​+sin2θ2​(cosθ1​cosθ2​+sinθ1​sinθ2​)+i(sinθ1​cosθ2​−cosθ1​sinθ2​)​]=r2​r1​​[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​


正题

某些奇奇怪怪的读法。

  1. f a ˉ   f a ˉ   t a ˇ f\bar{a}\ f\bar{a}\ t\check{a} faˉ faˉ taˇ
  2. F a s t   F a s t   T L E Fast\ Fast\ TLE Fast Fast TLE

Cooley-Tukey \texttt{Cooley-Tukey} Cooley-Tukey 算法

%   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−1​ai​xi  按照下标奇偶性分类,则有
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​+a2​x2+a4​x4+…+an−2​xn−2)+(a1​x+a3​x3+a5​x5+…+an−1​xn−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​+a2​x+a4​x2+…+an−2​x2n​−1A2​(x)=a1​+a3​x+a5​x2+…+an−1​x2n​−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​)+wnk​A2​(wn2k​)=A1​(w2n​k​)+wnk​A2​(w2n​k​)​

%   将 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+2n​​A2​(wn2k+n​)=A1​(wn2k​∗wnn​)−wnk​A2​(wn2k​∗wnn​)=A1​(wn2k​)−wnk​A2​(wn2k​)=A1​(w2n​k​)−wnk​A2​(w2n​k​)​  发现了什么?只相差一个符号!有什么作用?那么当我们在枚举第一个式子的时候,我们可以 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) Θ(大常数+nlog2​n)这个大常数是哪里来的呢?复数运算!复数的乘法是非常慢的。因此如果不是数据非常大,应尽量避免使用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​+a1​wn1​+a2​(wn1​)2+⋯+an​(wn1​)n=a0​+a1​wn1​+a2​wn2​+⋯+an​wnn​=j=0∑n−1​aj​(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−1​yi​(wn−k​)i=i=0∑n−1​(j=0∑n−1​aj​(wni​)j)(wn−k​)i=i=0∑n−1​(j=0∑n−1​aj​(wnj​)i)(wn−k​)i=i=0∑n−1​(j=0∑n−1​aj​(wnj​)i(wn−k​)i)=i=0∑n−1​j=0∑n−1​aj​(wnj​)i(wn−k​)i=i=0∑n−1​j=0∑n−1​aj​(wnj−k​)i=j=0∑n−1​aj​(i=0∑n−1​(wnj−k​)i)=j=0∑n−1​aj​(i=0∑n−1​(wnj−k​)i)​  根据性质5,有:

  1. 当 j ≠ k j\ne k j​=k 时,即 j − k ≠ 0 j-k\ne 0 j−k​=0,此时 ∑ i = 0 n − 1 ( w n j − k ) i = 0 \sum\limits_{i=0}^{n-1}(w_n^{j-k})^i=0 i=0∑n−1​(wnj−k​)i=0,共有 n − 1 n-1 n−1 种情况
  2. 当 j = k j=k j=k 时,即 j − k = 0 j-k=0 j−k=0,此时 ∑ i = 0 n − 1 ( w n j − k ) i = n \sum\limits_{i=0}^{n-1}(w_n^{j-k})^i=n i=0∑n−1​(wnj−k​)i=n,共有 1 1 1 种情况。
    ∴ c k = n × a k \therefore c_k=n\times a_k ∴ck​=n×ak​

%   因而有
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(不小的常数+Nlog2​N)
  下面是完整的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(Nlog2​N)的空间是开在栈空间里的,因此可能会爆栈!所以有了下一小节。

迭代 Cooley-Tukey \texttt{Cooley-Tukey} Cooley-Tukey 算法

%   观察一下原下标和新下标的二进制,发现了什么规律?
  规律不是很好阐述,这里定义 x x x 的 l l l 位翻转为 x x x 在长度为 log ⁡ 2 l − 1 \log_2l-1 log2​l−1 的二进制下的翻转1
  举个例子: 1 1 1 的 8 8 8 位翻转为 4 4 4。
   1 1 1 的长度为 log ⁡ 2 8 − 1 = 3 \log_28-1=3 log2​8−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(Nlog2​N)。
%   因此如果求翻转序列的复杂度比 O ( N l o g 2 N ) O(Nlog_2N) O(Nlog2​N) 还要高那就完了。

暴力

%   时间复杂度: O ( N l o g 2 N ) O(Nlog_2N) O(Nlog2​N)
%   空间复杂度: O ( n ) O(n) O(n)
%   代码复杂度:小
%   代码:下面的 l = log ⁡ 2 l i m i t l=\log_2 limit l=log2​limit,实现上可以在求 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

Luogu 1919

构造

[MUTC2013]idiots

BZOJ3513

容斥原理



  1. 其中 l = 2 x , x ∈ Z l=2^x,x\in\Z l=2x,x∈Z, log ⁡ \log log的优先级比减法高。 ↩︎