FFT
- 多项式
- 多项式乘法
- 复数及运算
- 导数
- 泰勒公式及展开式
- 欧拉公式
- 单位根
- FFT
- Code
- IFFT
多项式
我们从课本中可以知道,一个 n − 1 n-1 n−1 次的多项式可以写成 a 0 + a 1 x + a 2 x 2 + a 3 x 3 + ⋯ + a n − 1 x n − 1 a_{0}+a_{1}x+a_{2}x^2+a_{3}x^3+\dots+a_{n-1}x^{n-1} a0+a1x+a2x2+a3x3+⋯+an−1xn−1
用高级一点的表示法就是:
一个 n − 1 n-1 n−1 次 n n n 项多项式 f ( x ) f(x) f(x) 可以表示为:
f ( x ) = ∑ i = 0 n − 1 a i x i f(x)=\sum_{i=0}^{n-1}a_ix^i f(x)=∑i=0n−1aixi
所以我们可以用每一项的系数来表示 f ( x ) f(x) f(x),即
f ( x ) = { a 0 , a 1 , a 2 , … , a n − 1 } f(x)=\{a_0,a_1,a_2,\dots,a_{n-1}\} f(x)={a0,a1,a2,…,an−1}
俗称系数表示法。
但是,当我们把这个多项式看成一个函数时,我们将 n n n 个不同的 x x x 带入,会得到 n n n 个不一样的 y y y,而我们就得到了 n n n 个点的坐标,它们确定唯一的多项式 f ( x ) f(x) f(x)
证明很简单,直接待定系数法就可以解决。
那么 f ( x ) f(x) f(x) 还可以用这 n n n 个点表示:
f ( x ) = { ( x 0 , y 0 ) , ( x 1 , y 1 ) , ( x 2 , y 2 ) , … , ( x n − 1 , y n − 1 ) } f(x)=\{(x_0,y_0),(x_1,y_1),(x_2,y_2),\dots,(x_{n-1},y_{n-1})\} f(x)={(x0,y0),(x1,y1),(x2,y2),…,(xn−1,yn−1)}
(其中 y i = f ( x i ) y_i=f(x_i) yi=f(xi))
俗称点值表示法。
多项式乘法
对于两个系数表示法表示的多项式 f ( x ) f(x) f(x) 和 g ( x ) g(x) g(x),如果要计算它们的乘积 h ( x ) h(x) h(x),我们需要枚举 f ( x ) f(x) f(x) 的每一项和 g ( x ) g(x) g(x) 的每一项相乘,最后累计,时间复杂度 O ( n 2 ) O(n^2) O(n2)
那么换成点值表示呢?时间复杂度为 O ( n ) O(n) O(n)
Why?
设两个多项式 f ( x ) f(x) f(x) 和 g ( x ) g(x) g(x) 分别为:
f ( x ) = { ( x 0 , y 1 0 ) , ( x 1 , y 1 1 ) , ( x 2 , y 1 2 ) , … , ( x n − 1 , y 1 n − 1 ) } f(x)=\{(x_0,y1_0),(x_1,y1_1),(x_2,y1_2),\dots,(x_{n-1},y1_{n-1})\} f(x)={(x0,y10),(x1,y11),(x2,y12),…,(xn−1,y1n−1)}
(其中 y 1 i = f ( x i ) y1_i=f(x_i) y1i=f(xi))
g ( x ) = { ( x 0 , y 2 0 ) , ( x 1 , y 2 1 ) , ( x 2 , y 2 2 ) , … , ( x n − 1 , y 2 n − 1 ) } g(x)=\{(x_0,y2_0),(x_1,y2_1),(x_2,y2_2),\dots,(x_{n-1},y2_{n-1})\} g(x)={(x0,y20),(x1,y21),(x2,y22),…,(xn−1,y2n−1)}
(其中 y 2 i = g ( x i ) y2_i=g(x_i) y2i=g(xi))
那么结果 h ( x ) h(x) h(x) 为:
h ( x ) = { ( x 0 , y 2 0 ⋅ y 1 0 ) , ( x 1 , y 2 1 ⋅ y 1 1 ) , ( x 2 , y 2 2 ⋅ y 1 2 ) , … , ( x n − 1 , y 2 n − 1 ⋅ y 1 n − 1 ) } h(x)=\{(x_0,y2_0 \cdot y1_0),(x_1,y2_1 \cdot y1_1),(x_2,y2_2 \cdot y1_2),\dots,(x_{n-1},y2_{n-1} \cdot y1_{n-1})\} h(x)={(x0,y20⋅y10),(x1,y21⋅y11),(x2,y22⋅y12),…,(xn−1,y2n−1⋅y1n−1)}
好,我们就结束了
但是怎么将系数表示法转换为点值表示法呢?
复数及运算
高中课本里讲过复数,不过不知道也没关系,这里会讲。
定义:复数单位 i i i 满足 i 2 = − 1 i^2=-1 i2=−1,则形如 a + b i a+bi a+bi 的数为复数。 ( a , b ∈ R ) (a,b\in R) (a,b∈R)
其中 a a a 叫做实部, b b b 叫做虚部, i i i 为虚数单位。
所以有 − 7 = 7 i \sqrt{-7}=\sqrt7i −7=7i
类似平面直角坐标系,在复平面直角坐标系中,复数为其中的一个点。
如下图:
点(3,2)表示的复数为 3 + 2 i 3+2i 3+2i
它也可以表示为( 1 3 , θ \sqrt 13,\theta 13,θ)
定义:一个复数的模为它到原点的距离。
及复数 z = a + b i z=a+bi z=a+bi 的模记为 ∣ z ∣ = a 2 + b 2 |z|=\sqrt{a^2+b^2} ∣z∣=a2+b2
一个复数的共轭复数为其虚部取反。
及复数 z = a + b i z=a+bi z=a+bi 的共轭复数为 z ˉ = a − b i \bar{z}=a-bi zˉ=a−bi, z z ˉ = a 2 + b 2 z \bar{z} = a^2 +b^2 zzˉ=a2+b2
其实就跟共轭根式差不多。
复数的运算其实跟实数的运算差不多。
设 z 1 = a + b i z_1 = a + bi z1=a+bi, z 2 = c + d i z_2 = c + di z2=c+di
加减法:实部和虚部分别相加减:
z 1 ± z 2 = ( a ± c ) + ( b ± d ) i z_1 \pm z_2 = (a \pm c) + (b \pm d)i z1±z2=(a±c)+(b±d)i
乘法:硬算出奇迹(拆完后合并同类项)
z 1 z 2 = ( a c − b d ) + ( b c + a d ) i z_1z_2 = (ac - bd) + (bc + ad)i z1z2=(ac−bd)+(bc+ad)i
( a 1 , θ 1 ) ( a 2 , θ 2 ) = ( a 1 a 2 , θ 1 + θ 2 ) (a_1,\theta_1)(a_2,\theta_2)=(a_1a_2,\theta_1+\theta_2) (a1,θ1)(a2,θ2)=(a1a2,θ1+θ2)
除法:化简
z 1 z 2 = a + b i c + d i = ( a + b i ) ( c − d i ) ( c + d i ) ( c − d i ) = a c − b d i 2 + b c i − a d i c 2 − d 2 i 2 = ( a c + b d ) + ( b c − a d ) i c 2 + d 2 \dfrac{z_1}{z_2} = \dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{ac-bdi^2+bci-adi}{c^2-d^2i^2}=\dfrac{(ac+bd)+(bc-ad)i}{c^2+d^2} z2z1=c+dia+bi=(c+di)(c−di)(a+bi)(c−di)=c2−d2i2ac−bdi2+bci−adi=c2+d2(ac+bd)+(bc−ad)i
导数
和初中的 y = … y=\dots y=… 不一样,我们用 f ( x ) f(x) f(x) 来表示一个关于 x x x 的函数。
一般地,已知函数 y = f ( x ) y=f(x) y=f(x), x 0 , x 1 x_0,x_1 x0,x1 是其定义域内不同两点,记 Δ x = x 1 − x 0 , Δ y = y 1 − y 0 = f ( x 1 ) − f ( x 0 ) = f ( x 0 + Δ x ) − f ( x 0 ) \Delta x=x_1-x_0,\Delta y=y_1-y_0=f(x_1)-f(x_0)=f(x_0+\Delta x)-f(x_0) Δx=x1−x0,Δy=y1−y0=f(x1)−f(x0)=f(x0+Δx)−f(x0),则当 Δ x ≠ 0 \Delta x\neq 0 Δx=0 时,商 f ( x 0 + Δ x ) − f ( x 0 ) Δ x = Δ y Δ x \frac{f(x_0+\Delta x)-f(x_0)}{\Delta x}=\frac{\Delta y}{\Delta x} Δxf(x0+Δx)−f(x0)=ΔxΔy 称作函数 y = f ( x ) y=f(x) y=f(x)在区间 $[x_0,x_0 +_\Delta x] $ 上的平均变化率。
当 Δ x \Delta x Δx 趋近于 0 0 0 时,平均变化率 Δ y Δ x = f ( x 0 + Δ x ) − f ( x 0 ) Δ x \frac{\Delta y}{\Delta x}=\frac{f(x_0+\Delta x)-f(x_0)}{\Delta x} ΔxΔy=Δxf(x0+Δx)−f(x0) 趋近于一个常数 l l l,那么常数 l l l 称为函数 f ( x ) f(x) f(x) 在点 x 0 x_0 x0 的瞬时变化率。
趋近于可以用符号 → \to → 表示,所以上面那句话可以这样表示:
当 Δ x → 0 \Delta x \to 0 Δx→0 时, f ( x 0 + Δ x ) − f ( x 0 ) Δ x → l \frac{f(x_0+\Delta x)-f(x_0)}{\Delta x} \to l Δxf(x0+Δx)−f(x0)→l
或记作:
lim Δ x → 0 f ( x 0 + Δ x ) − f ( x 0 ) Δ x = l \lim_{\Delta x \to 0} \frac{f(x_0+\Delta x)-f(x_0)}{\Delta x}=l limΔx→0Δxf(x0+Δx)−f(x0)=l
函数在 x 0 x_0 x0 的瞬时变化率,通常称为 f ( x ) f(x) f(x) 在 x = x 0 x=x_0 x=x0 处的导数,并记作 f ′ ( x 0 ) f'(x_0) f′(x0)。
注意,是 f ′ f' f′。
那么对于函数 f ( x ) f(x) f(x) 怎么求导呢?
如下:
{ ( x a ) ′ = a x a − 1 ( c ) ′ = 0 ( e x ) ′ = e x ( ln x ) ′ = 1 x ( sin x ) ′ = cos x ( cos x ) ′ = sin x } \begin{Bmatrix} (x^a)'=ax^{a-1} \\ (c)'=0 \\ (e^x)'=e^x \\ (\ln x)'=\frac{1}{x} \\ (\sin x)'=\cos x \\ (\cos x)'=\sin x \end{Bmatrix} ⎩ ⎨ ⎧(xa)′=axa−1(c)′=0(ex)′=ex(lnx)′=x1(sinx)′=cosx(cosx)′=sinx⎭ ⎬ ⎫
例如 f ( x ) = x 2 f(x)=x^2 f(x)=x2,则 f ′ ( x ) = 2 x f'(x)=2x f′(x)=2x
最后再说一点 f ′ ( x ) f'(x) f′(x) 表示函数 f ( x ) f(x) f(x) 在 x x x 取值为几时切线的斜率。
泰勒公式及展开式
泰勒公式:
f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) 2 ! ( x − x 0 ) 2 + ⋯ + f ( n ) ( x 0 ) n ! ( x − x 0 ) n + R n + 1 ( x ) f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\dotsb +\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_{n+1}(x) f(x)=f(x0)+f′(x0)(x−x0)+2!f′′(x0)(x−x0)2+⋯+n!f(n)(x0)(x−x0)n+Rn+1(x)
其中 R R R 为误差。
那么当 x 0 = 0 x_0=0 x0=0 时,得到麦克劳林公式:
f ( x ) = f ( 0 ) + f ′ ( 0 ) x + f ′ ′ ( 0 ) 2 ! x 2 + ⋯ + f ( n ) ( 0 ) n ! x n + R n + 1 ( x ) f(x)=f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+\dotsb +\frac{f^{(n)}(0)}{n!}x^n+R_{n+1}(x) f(x)=f(0)+f′(0)x+2!f′′(0)x2+⋯+n!f(n)(0)xn+Rn+1(x)
所以我们就有了以下几个展开式:
{ e x = 1 + x + x 2 2 ! + x 3 3 ! + ⋯ + x n n ! + ⋯ ln ( 1 + x ) = x − x 2 2 + x 3 3 + ⋯ + ( − 1 ) n − 1 x n n + ⋯ sin x = ∑ k = 0 ∞ ( − 1 ) k x 2 k + 1 ( 2 k + 1 ) ! cos x = ∑ k = 0 ∞ ( − 1 ) k x 2 k ( 2 k ) ! } \begin{Bmatrix} e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\dotsb+\frac{x^n}{n!}+\dotsb \\ \ln(1+x)=x-\frac{x^2}{2}+\frac{x^3}{3}+\dotsb+(-1)^{n-1}\frac{x^n}{n}+\dotsb \\ \sin x=\sum_{k=0}^{\infty}(-1)^k\frac{x^{2k+1}}{(2k+1)!} \\ \cos x=\sum_{k=0}^{\infty}(-1)^k\frac{x^{2k}}{(2k)!} \end{Bmatrix} ⎩ ⎨ ⎧ex=1+x+2!x2+3!x3+⋯+n!xn+⋯ln(1+x)=x−2x2+3x3+⋯+(−1)n−1nxn+⋯sinx=∑k=0∞(−1)k(2k+1)!x2k+1cosx=∑k=0∞(−1)k(2k)!x2k⎭ ⎬ ⎫
欧拉公式
接上文,展开 e i x e^{ix} eix,得:
e i x = 1 + i x − x 2 2 ! − i x 3 3 ! + x 4 4 ! + i x 5 5 ! − x 6 6 ! − i x 7 7 ! + ⋯ ) e^{ix}=1+ix-\frac{x^2}{2!}-i\frac{x^3}{3!}+\frac{x^4}{4!}+i\frac{x^5}{5!}-\frac{x^6}{6!}-i\frac{x^7}{7!}+\dotsb) eix=1+ix−2!x2−i3!x3+4!x4+i5!x5−6!x6−i7!x7+⋯)
e i x = ( 1 − x 2 2 ! + x 4 4 ! − x 6 6 ! + ⋯ ) + i ( x − x 3 3 ! + x 5 5 ! − 7 3 7 ! + ⋯ ) e^{ix}=(1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\dotsb) +i(x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{7^3}{7!}+\dotsb) eix=(1−2!x2+4!x4−6!x6+⋯)+i(x−3!x3+5!x5−7!73+⋯)
e i x = cos x + i sin x e^{ix}=\cos x+i\sin x eix=cosx+isinx
将 x = π x=\pi x=π 带入:
e i π = cos π + i sin π = − 1 e^{i\pi}=\cos \pi + i\sin \pi=-1 eiπ=cosπ+isinπ=−1
所以欧拉恒等式出来:
e i π + 1 = 0 e^{i\pi}+1=0 eiπ+1=0
单位根
设 n n n 为正整数,若 x n = 1 x^n=1 xn=1 ,则 x x x 是 n n n 次单位根。
下文默认 n = 2 m , m n=2^m,m n=2m,m为正整数。
根据上文的欧拉公式,我们可知 e 2 k π i = cos 2 k π + i sin 2 k π = 1 e^{2k \pi i}=\cos 2k\pi+i\sin 2k\pi=1 e2kπi=cos2kπ+isin2kπ=1,故 n n n 次单位根分别是: e 2 k π i n = cos 2 k π n + i sin 2 k π n , ( k = 0 , 1 , 2 , … , n − 1 ) e^{\frac{2k\pi i}{n}}=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n},(k=0,1,2,\dots,n-1) en2kπi=cosn2kπ+isinn2kπ,(k=0,1,2,…,n−1)。
接上文复数的知识,每个单位根在复平面上的坐标为 ( cos 2 k π n , sin 2 k π n ) (\cos \frac{2k\pi}{n},\sin \frac{2k\pi}{n}) (cosn2kπ,sinn2kπ)。
如图:
在一个标准复平面坐标系上,以 1 1 1 为半径画一个圆, n n n次单位根可以想象为把这个圆平分为 n n n 分中每一个点。
我们用 ω n 1 , ω n 2 , ω n 3 , … , ω n n \omega_n^1,\omega_n^2,\omega_n^3,\dots,\omega_n^n ωn1,ωn2,ωn3,…,ωnn 来表示 n n n 次单位根的每一个根,及方程的每一个解。
所以有 ω n k = cos k 2 π n + i sin k 2 π n \omega_n^k=\cos k\frac{2\pi}{n}+i\sin k\frac{2\pi}{n} ωnk=coskn2π+isinkn2π
如当 n = 4 n=4 n=4时, ω 4 = 1 , − 1 , i , − 1 \omega_4=1,-1,i,-1 ω4=1,−1,i,−1
如当 n = 1 n=1 n=1 时, ω 1 = 1 \omega_1=1 ω1=1
单位根的性质:
ω 2 n 2 k = ω n k = − ω n k + n 2 \omega_{2n}^{2k}=\omega_{n}^{k}=-\omega _n^{k+\frac{n}{2}} ω2n2k=ωnk=−ωnk+2n
ω n 0 = ω n n = 1 \omega_n^0=\omega_n^n=1 ωn0=ωnn=1
ω 2 n 2 k = cos 2 k 2 π 2 n + i sin 2 k 2 π 2 n = ω n k \omega_{2n}^{2k}=\cos 2k \frac{2\pi}{2n}+i\sin 2k \frac{2\pi}{2n}=\omega_n^k ω2n2k=cos2k2n2π+isin2k2n2π=ωnk
ω n n 2 = cos n 2 ⋅ 2 π n + i sin n 2 ⋅ 2 π n = cos π + i sin π = − 1 \omega_n^{\frac{n}{2}}=\cos \frac{n}{2}\cdot \frac{2\pi}{n}+i\sin \frac{n}{2}\cdot \frac{2\pi}{n}=\cos \pi+i\sin \pi=-1 ωn2n=cos2n⋅n2π+isin2n⋅n2π=cosπ+isinπ=−1
FFT
设 n n n 为偶数, k < n 2 k<\frac{n}{2} k<2n。
设
A ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + ⋯ + a n − 1 x n − 1 A(x)=a_{0}+a_{1}x+a_{2}x^2+a_{3}x^3+\dots+a_{n-1}x^{n-1} A(x)=a0+a1x+a2x2+a3x3+⋯+an−1xn−1
A 1 ( x ) = a 0 + a 2 x + a 4 x 2 + ⋯ + a n − 2 x n 2 − 1 A_1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1} A1(x)=a0+a2x+a4x2+⋯+an−2x2n−1
A 2 ( x ) = a 1 + a 3 x + a 5 x 2 + ⋯ + a n − 1 x n 2 − 1 A_2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1} A2(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)
A ( ω n k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k}) A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)
A ( ω n k + n 2 ) = A 1 ( ω n 2 k + n ) + ω n k + n 2 A 2 ( ω n 2 k + n ) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k ) A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k}) A(ωnk+2n)=A1(ωn2k+n)+ωnk+2nA2(ωn2k+n)=A1(ωn2k)−ωnkA2(ωn2k)
发现了什么?没错,只有后面一坨是相反的,这说明我们可以 O ( 1 ) O(1) O(1) 通过第一个式子得出第二个式子
所以 k k k 只要枚举前 n 2 \frac{n}{2} 2n 个数就可以了,问题缩小了一半,而缩小后还能继续缩小,所以我们用类似线段树的分治(递归)来解决它,时间复杂度 O ( n log n ) O(n\log n) O(nlogn)
Code
void fft(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];fft(limit>>1,a1,type);fft(limit>>1,a2,type);complex Wn=complex(cos(2.0*Pi/limit),type*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];
}
IFFT
但是光有点值还是不够的,我们还要转换为数值,这就是 IFFT。
至于还有一个迭代优化。
我们下篇博客再讲。