【算法】FFT-1(递归实现)(不包括IFFT)

多项式

我们从课本中可以知道,一个 n − 1 n-1 n1 次的多项式可以写成 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++an1xn1

用高级一点的表示法就是:
一个 n − 1 n-1 n1 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=0n1aixi

所以我们可以用每一项的系数来表示 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,,an1}

俗称系数表示法。

但是,当我们把这个多项式看成一个函数时,我们将 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),,(xn1,yn1)}

(其中 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),,(xn1,y1n1)}

(其中 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),,(xn1,y2n1)}

(其中 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,y20y10),(x1,y21y11),(x2,y22y12),,(xn1,y2n1y1n1)}

好,我们就结束了

但是怎么将系数表示法转换为点值表示法呢?


复数及运算

高中课本里讲过复数,不过不知道也没关系,这里会讲。

定义:复数单位 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,bR)

其中 a a a 叫做实部, b b b 叫做虚部, i i i 为虚数单位。

所以有 − 7 = 7 i \sqrt{-7}=\sqrt7i 7 =7 i

类似平面直角坐标系,在复平面直角坐标系中,复数为其中的一个点。

如下图:

点(3,2)表示的复数为 3 + 2 i 3+2i 3+2i

它也可以表示为( 1 3 , θ \sqrt 13,\theta 1 3,θ

定义:一个复数的模为它到原点的距离。

及复数 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ˉ=abi 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=(acbd)+(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)(cdi)(a+bi)(cdi)=c2d2i2acbdi2+bciadi=c2+d2(ac+bd)+(bcad)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=x1x0,Δy=y1y0=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 Δx0 时, 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Δx0Δ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)=axa1(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)(xx0)+2!f′′(x0)(xx0)2++n!f(n)(x0)(xx0)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)=x2x2+3x3++(1)n1nxn+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+ix2!x2i3!x3+4!x4+i5!x56!x6i7!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=(12!x2+4!x46!x6+)+i(x3!x3+5!x57!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 e=cosπ+isinπ=1

所以欧拉恒等式出来:

e i π + 1 = 0 e^{i\pi}+1=0 e+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=cos2+isin2=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=cosn2+isinn2,(k=0,1,2,,n1)

接上文复数的知识,每个单位根在复平面上的坐标为 ( cos ⁡ 2 k π n , sin ⁡ 2 k π n ) (\cos \frac{2k\pi}{n},\sin \frac{2k\pi}{n}) (cosn2,sinn2)
如图:

在一个标准复平面坐标系上,以 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=cos2nn2π+isin2nn2π=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++an1xn1
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++an2x2n1
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++an1x2n1
所以有:
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。
至于还有一个迭代优化。
我们下篇博客再讲。

  • 62
    点赞
  • 31
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值