超详细易懂FFT(快速傅里叶变换)及代码实现_傅立叶变换编程-程序员宅基地

技术标签: 数论  

前言

昨天学了一晚上,终于搞懂了FFT。希望能写一篇清楚易懂的题解分享给大家,也进一步加深自己的理解。
FFT算是数论中比较重要的东西,听起来就很高深的亚子。但其实学会了(哪怕并不能完全理解),会实现代码,并知道怎么灵活运用 (背板子) 就行。接下来进入正题。

定义

FFT(Fast Fourier Transformation),中文名快速傅里叶变换,是离散傅氏变换的快速算法,它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。
而在信奥中,一般用来加速多项式乘法
朴素高精度乘法的时间为 O ( n 2 ) O(n^2) O(n2),但FFT能将时间复杂度降到 O ( n l o g 2 n ) O(nlog_2n) O(nlog2n)
学习FFT之前,需要了解一些有关复数和多项式的知识。

有关知识

多项式的两种表示方法

系数表示法

F [ x ] = y = a 0 x 0 + a 1 x 1 + a 2 x 2 + . . . . . . a n x n F[x]=y=a_0x^0+a_1x^1+a_2x^2+......a_nx^n F[x]=y=a0x0+a1x1+a2x2+......anxn
{ a 0 , a 1 , a 2 , . . . , a n a_0,a_1,a_2,...,a_n a0,a1,a2,...,an} 是这个多项式每一项的系数,所以这是多项式的系数表示法

点值表示法

在函数图像中, F [ x ] F[x] F[x]这个多项式可以被n个点唯一确定,即代入n个点作为 x x x,分别解出对应的 y y y,得到n个式子。把这n条式子联立起来成为一个有n条方程的n元方程组,每一项的系数都可以解出来.(可类比二元一次方程)
也就是说,使用{ ( x 0 , f [ x 0 ] ) (x_0,f[x_0]) (x0,f[x0]), ( x 1 , f [ x 1 ] ) (x_1,f[x_1]) (x1,f[x1]),…, ( x n , f [ x n ] ) (x_n,f[x_n]) (xn,f[xn])}就可以完整描述出这个多项式,这就是 多项式的点值表示法

多项式相乘

设两个多项式分别为 f ( x ) f(x) f(x), g ( x ) g(x) g(x),我们要把这两个多项式相乘 (即求卷积)。
如果用系数表示法:
我们要枚举 f f f的每一位的系数与 g g g的每一位的系数相乘,多项式乘法时间复杂度 O ( n 2 ) O(n^2) O(n2),这也是我们所熟知的高精度乘法的原理。
如果用点值表示法:
f [ x ] f[x] f[x]={ ( x 0 , f [ x 0 ] ) (x_0,f[x_0]) (x0,f[x0]), ( x 1 , f [ x 1 ] ) (x_1,f[x_1]) (x1,f[x1]),…, ( x n , f [ x n ] ) (x_n,f[x_n]) (xn,f[xn])}
g [ x ] g[x] g[x]={ ( x 0 , g [ x 0 ] ) (x_0,g[x_0]) (x0,g[x0]), ( x 1 , g [ x 1 ] ) (x_1,g[x_1]) (x1,g[x1]),…, ( x n , g [ x n ] ) (x_n,g[x_n]) (xn,g[xn])}
f [ x ] ∗ g [ x ] f[x]*g[x] f[x]g[x]={ ( x 0 , f [ x 0 ] ∗ g [ x 0 ] ) (x_0,f[x_0]*g[x_0]) (x0,f[x0]g[x0]), ( x 1 , f [ x 1 ] ∗ g [ x 1 ] ) (x_1,f[x_1]*g[x_1]) (x1,f[x1]g[x1]),…, ( x n , f [ x n ] ∗ g [ x n ] ) (x_n,f[x_n]*g[x_n]) (xn,f[xn]g[xn])}
我们可以发现,如果两个多项式取相同的 x x x,得到不同的 y y y值,那么只需要 y y y值对应相乘就可以了!
复杂度只有枚举 x x x O ( n ) O(n) O(n)
那么问题转换为将多项式系数表示法转化成点值表示法。
朴素系数转点值的算法叫DFT(离散傅里叶变换),优化后为FFT(快速傅里叶变换)点值转系数的算法叫IDFT(离散傅里叶逆变换),优化后为IFFT(快速傅里叶逆变换)。之后我会分别介绍。

卷积

其实不理解卷积也没关系,但这里顺便提一下,可以跳过的
卷积与傅里叶变换有着密切的关系。利用一点性质,即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换,能使傅里叶分析中许多问题的处理得到简化。

F ( g ( x ) ∗ f ( x ) ) = F ( g ( x ) ) F ( f ( x ) ) F(g(x)*f(x)) = F(g(x))F(f(x)) F(g(x)f(x))=F(g(x))F(f(x))

其中F表示的是傅里叶变换

复数

高中数学会详细讲解,知道的可以跳过这一部分,没学过也没关系,看以下内容应该能很清楚的理解。

1.定义

数集拓展到实数范围内,仍有些运算无法进行。比如判别式小于0的一元二次方程仍无解,因此将数集再次扩充,达到复数范围。
复数 z z z被定义为二元有序实数对 ( a , b ) (a,b) (a,b),记为 z = a + b i z=a+bi z=a+bi,这里 a a a b b b是实数,规定 i i i是虚数单位。 ( i 2 = − 1 i^2=-1 i2=1 i = − 1 i=\sqrt{-1} i=1 )
对于复数 z = a + b i z=a+bi z=a+bi。实数 a a a称为复数z的实部(real part),记作 r e z = a rez=a rez=a.实数 b b b称为复数z的虚部(imaginary part)记作 Imz=b.
当虚部等于零时,这个复数可以视为实数
即当 b = 0 b=0 b=0时, z = a z=a z=a,这时复数成为实数;当且仅当 a = b = 0 a=b=0 a=b=0时,它是实数0;
当z的虚部不等于零时,实部等于零时,常称z为纯虚数
即当 a = 0 a=0 a=0 b ≠ 0 b≠0 b=0时, z = b i z=bi z=bi,我们就将其称为纯虚数。
将复数的实部与虚部的平方和的正的平方根的值称为该复数的模,记作 ∣ z ∣ ∣z∣ z
即对于复数z=a+bi,它的模为 ∣ z ∣ = ( a 2 + b 2 ) ∣z∣=\sqrt{(a^2+b^2)} z=(a2+b2)

2.复数的几何意义

直接两张图搞定√ (应该可以一目了然)
在这里插入图片描述
在这里插入图片描述

3.运算法则

加法法则: ( a + b i ) + ( c + d i ) = ( a + c ) + ( b + d ) i ; (a+bi)+(c+di)=(a+c)+(b+d)i; (a+bi)+(c+di)=(a+c)+(b+d)i;
减法法则: ( a + b i ) − ( c + d i ) = ( a − c ) + ( b − d ) i ; (a+bi)-(c+di)=(a-c)+(b-d)i; (a+bi)(c+di)=(ac)+(bd)i;
注:复数加减满足平行四边形法则

在这里插入图片描述
乘法法则: ( a + b i ) ⋅ ( c + d i ) = ( a c − b d ) + ( b c + a d ) i (a+bi)·(c+di)=(ac-bd)+(bc+ad)i (a+bi)(c+di)=(acbd)+(bc+ad)i
复数相乘一个重要法则:模长相乘,幅角相加。(这个定理很重要)
模长:这个向量的模长,即这个点到原点的距离。(不懂的可再看下向量的几何意义)。
幅角: 从原点出发、指向x轴正半轴的射线绕原点逆时针旋转至过这个点所经过的角。
在极坐标(可看成平面直角坐标系)下,复数可用模长r与幅角θ表示为 ( r , θ ) (r,θ) (r,θ)。对于复数 a + b i a+bi a+bi, r = ( a ² + b ² ) r=\sqrt{(a²+b²)} r=(a²+b²) θ = a r c t a n ( b / a ) θ=arctan(b/a) θ=arctan(b/a)。此时,复数相乘表现为模长相乘,幅角相加。
除法法则 ( a + b i ) ÷ ( c + d i ) = [ ( a c + b d ) / ( c ² + d ² ) ] + [ ( b c − a d ) / ( c ² + d ² ) ] i (a+bi)÷(c+di)=[(ac+bd)/(c²+d²)]+[(bc-ad)/(c²+d²)]i (a+bi)÷(c+di)=[(ac+bd)/(c²+d²)]+[(bcad)/(c²+d²)]i

4. 共轭复数

一个复数 z = a + b i z=a+bi z=a+bi共轭复数 a − b i a−bi abi(实部不变,虚部取反),记为 z ‾ = a − b i \overline{z}=a-bi z=abi
当复数模为1时(即|z|=1),与共轭复数互为倒数
证明: z ∗ z ‾ = a 2 − b 2 ∗ i 2 = a 2 + b 2 = ∣ z ∣ 2 = 1 z*\overline{z}=a^2-b^2*i^2=a^2+b^2=|z|^2=1 zz=a2b2i2=a2+b2=z2=1

FFT加速多项式乘法

由于多项式乘法用点值表示比用系数表示快的多,所以我们先要将系数表示法转化成点值表示法相乘,再将结果的点值表示法转化为系数表示法的过程。
第一个过程叫做FFT(快速傅里叶变换),第二个过程叫IFFT(快速傅里叶逆变换)
在讲这两个过程之前,首先了解一个概念:

单位根

复数 ω \omega ω满足 ω n = 1 \omega^n=1 ωn=1,称 ω \omega ω是n次单位根

怎么找单位根?

单位圆:圆心为原点、1为半径的圆
单位圆n等分,取这n个点(或点表示的向量)所表示的复数(即分别以这n个点的横坐标为实部、纵坐标为虚部,所构成的虚数),即为n次单位根
下图包含了当n=8时,所有的8次单位根,分别记为 ω 8 1 , ω 8 2 . . . . . , ω 8 8 \omega_8^1,\omega_8^2.....,\omega_8^8 ω81,ω82.....,ω88
(图中圆的半径是1,w表示 ω \omega ω,且下标8已省略)
图是我自己画的,可能有点丑QWQ
单位根图像byTrilarflagz
由此我们知道如何找单位根啦
从点(1,0)开始(即 ω n 1 \omega_n^1 ωn1),逆时针将这n个点从0开始编号,第k个点对应的虚数记作 ω n k \omega_n^k ωnk
由复数相乘法则:模长相乘幅角相加​ 可得:
( ω n 1 ) k = ω n k (\omega_n^1)^k=\omega_n^k (ωn1)k=ωnk
根据每个复数的幅角,可以计算出所对应的点/向量。 ω n k \omega_n^k ωnk 对应的点/向量是 ( c o s ⁡ k n 2 π , s i n ⁡ k n 2 π ) (cos⁡\frac kn2π,sin⁡\frac kn2π) (cosnk2π,sinnk2π),即为复数 c o s ⁡ k n 2 π + i ∗ s i n ⁡ k n 2 π cos⁡\frac kn2π+i *sin⁡\frac kn2π cosnk2π+isinnk2π

单位根的性质

建议记住,因为对之后的分析很重要!!

1. ω n k = ω 2 n 2 k \omega_n^k=\omega_{2n}^{2k} ωnk=ω2n2k
2. ω n k = − ω n k + n 2 \omega_n^k=-\omega_{n}^{k+\frac n 2} ωnk=ωnk+2n
3. ω n 0 = ω n n = 1 \omega_n^0=\omega_{n}^n=1 ωn0=ωnn=1

至于怎么证明,就是复数相乘时模长相乘幅角相加的原则。或者你直接观察图也可以很显然的得出结论。​

DFT(离散傅里叶变换)

对于任意多项式系数表示转点值表示,例如 F [ x ] = y = a 0 x 0 + a 1 x 1 + a 2 x 2 + . . . . . . + a n x n F[x]=y=a_0x^0+a_1x^1+a_2x^2+......+a_nx^n F[x]=y=a0x0+a1x1+a2x2+......+anxn ,可以随便取任意n个 x x x值代入计算,但这样时间复杂度是 O ( n 2 ) O(n^2) O(n2)
所以伟大数学家傅里叶取了一些特殊的点代入,从而进行优化。
他规定了点值表示中的 n n n x x x n n n个模长为1的复数。这 n n n个复数不是随机的,而是单位根
把上述的n个复数(单位根) ω n 0 , ω n 1 . . . . . , ω n n − 1 \omega_n^0,\omega_n^1.....,\omega_n^{n-1} ωn0,ωn1.....,ωnn1代入多项式,能得到一种特殊的点值表示,这种点值表示就叫DFT(离散傅里叶变换)

FFT(快速傅里叶变换)

虽然DFT能把多项式转换成点值但它仍然是暴力代入 n n n个数,复杂度仍然是O(n2),所以它只是快速傅里叶变换的朴素版。
所以我们要考虑利用单位根的性质,加速我们的运算,得到FFT(快速傅里叶变换)
对于多项式 A ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 A(x)=a_0+a_1x+a_2x^2+...+a_{n−1}x^{n−1} A(x)=a0+a1x+a2x2+...+an1xn1
将A(x)的每一项按照下标的奇偶分成两部分:
A ( x ) = a 0 + a 2 x 2 + . . . + a n − 2 x n − 2 + x ∗ ( a 1 + a 3 x 2 + . . . + a n − 1 x n − 2 ) A(x)=a_0+a_2x^2+...+a_{n−2}x^{n−2}+x*(a_1+a_3x^2+...+a_{n−1}x^{n−2}) A(x)=a0+a2x2+...+an2xn2+x(a1+a3x2+...+an1xn2)
设两个多项式 A 0 ( x ) A_0(x) A0(x) A 1 ( x ) A_1(x) A1(x),令:
A 0 ( x ) = a 0 x 0 + a 2 x 1 + . . . + a n − 2 x n / 2 − 1 A_0(x)=a_0x^0+a_2x^1+...+a_{n-2}x^{n/2-1} A0(x)=a0x0+a2x1+...+an2xn/21
A 1 ( x ) = a 1 x 0 + a 3 x 1 + . . . + a n − 1 x n / 2 − 1 A_1(x)=a_1x^0+a_3x^1+...+a_{n-1}x^{n/2-1} A1(x)=a1x0+a3x1+...+an1xn/21
显然, A ( x ) = A 0 ( x 2 ) + x ∗ A 1 ( x 2 ) A(x)=A_0(x^2)+x*A_1(x^2) A(x)=A0(x2)+xA1(x2)
假设 k < n k<n k<n,代入 x = ω n k x=ω_n^k x=ωnk(n次单位根):
A ( ω n k ) A(\omega_n^k) A(ωnk) = A 0 ( ω n 2 k ) + ω n k ∗ A 1 ( ω n 2 k ) =A_0(\omega_n^{2k})+\omega_n^{k}*A_1(\omega_n^{2k}) =A0(ωn2k)+ωnkA1(ωn2k)
= A 0 ( ω n 2 k ) + ω n k ∗ A 1 ( ω n 2 k ) =A_0(\omega_\frac n2^{k})+\omega_n^{k}*A_1(\omega_\frac n 2^{k}) =A0(ω2nk)+ωnkA1(ω2nk)

A ( ω n k + n 2 ) = A 0 ( ω n 2 k + n ) + ω n k + n 2 ∗ A 1 ( ω n 2 k + n ) A(\omega_n^{k+\frac n 2})=A_0(\omega_n^{2k+n})+\omega_n^{k+\frac n 2}*A_1(\omega_n^{2k+n}) A(ωnk+2n)=A0(ωn2k+n)+ωnk+2nA1(ωn2k+n)
= A 0 ( ω n 2 k ) − ω n k ∗ A 1 ( ω n 2 k ) =A_0(\omega_\frac n2^{k})-\omega_n^{k}*A_1(\omega_\frac n 2^{k}) =A0(ω2nk)ωnkA1(ω2nk)

考虑A1(x)和A2(x)分别在 ( ω n 2 1 , ω n 2 2 , ω n 2 3 , . . . , ω n 2 n 2 − 1 ) (\omega_\frac n 2^{1},\omega_\frac n 2^{2},\omega_\frac n 2^{3},...,\omega_\frac n 2^{\frac n 2-1}) (ω2n1,ω2n2,ω2n3,...,ω2n2n1)的点值表示已经求出,就可以O(n)求出A(x)在 ( ω n 1 , ω n 2 , ω n 3 , . . . , ω n n − 1 ) (\omega_n ^{1},\omega_n ^{2},\omega_n ^{3},...,\omega_n ^{n-1}) (ωn1,ωn2,ωn3,...,ωnn1)处的点值表示。这个操作叫蝴蝶变换
而A1(x)和A2(x)是规模缩小了一半的子问题,所以不断向下递归分治。当n=1的时候返回。
:这个过程一定要求每层都可以分成两大小相等的部分,所以多项式最高次项一定是2的幂,不是的话直接在最高次项补零QAQ。
时间复杂度 O ( n l o g 2 n ) O(nlog_2n) O(nlog2n)

IFFT(快速傅里叶逆变换)

我们已经将两个多项式从系数表示法转化成点值表示法相乘后,还要将结果从点值表示法转化为系数表示法,也就是IFFT(快速傅里叶逆变换)
首先思考一个问题,为什么要把 ω n k \omega_n^k ωnk(单位根)作为x代入?
当然是因为离散傅里叶变换特殊的性质,而这也和IFFT有关。
一个重要结论
把多项式A(x)的离散傅里叶变换结果作为另一个多项式B(x)的系数,取单位根的倒数即 ω n 0 , ω n − 1 . . . . . , ω n 1 − n \omega_n^0,\omega_n^{-1}.....,\omega_n^{1-n} ωn0,ωn1.....,ωn1n作为x代入B(x),得到的每个数再除以n,得到的是A(x)的各项系数,这就实现了傅里叶变换的逆变换了。相当于在FFT基础上再搞一次FFT。
证明(个人觉得写的非常清楚,不想看的跳过吧)~~

( y 0 , y 1 , y 2 , . . . , y n − 1 ) (y_0,y_1,y_2,...,y_{n−1}) (y0,y1,y2,...,yn1)为多项式
A ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 A(x)=a_0+a_1x+a_2x^2+...+a_{n−1}x^{n−1} A(x)=a0+a1x+a2x2+...+an1xn1的离散傅里叶变换。
设多项式 B ( x ) = y 0 + y 1 x + y 2 x 2 + . . . + y n − 1 x n − 1 B(x)=y_0+y_1x+y_2x^2+...+y_{n−1}x^{n−1} B(x)=y0+y1x+y2x2+...+yn1xn1
把离散傅里叶变换的 ω n 0 , ω n 1 . . . . . , ω n n − 1 \omega_n^0,\omega_n^1.....,\omega_n^{n-1} ωn0,ωn1.....,ωnn1这n个单位根的倒数,即 ω n 0 , ω n − 1 . . . . . , ω n 1 − n \omega_n^0,\omega_n^{-1}.....,\omega_n^{1-n} ωn0,ωn1.....,ωn1n作为x代入 B ( x ) B(x) B(x), 得到一个新的离散傅里叶变换 ( z 0 , z 1 , z 2 , . . . , z n − 1 ) (z_0,z_1,z_2,...,z{n−1}) (z0,z1,z2,...,zn1)
z k z_k zk= ∑ i = 0 n − 1 y i ( ω n − k ) i \sum_{i=0}^{n−1}y_i(ω_n^{-k})^i i=0n1yi(ωnk)i
= ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ∗ ( ω n i ) j ) ( ω n − k ) i \sum_{i=0}^{n−1}(\sum_{j=0}^{n-1}a_j*(\omega_n^i)^j)(ω_n^{-k})^i i=0n1(j=0n1aj(ωni)j)(ωnk)i
= ∑ j = 0 n − 1 a j ∗ ( ∑ i = 0 n − 1 ( ω n i ) j − k ) \sum_{j=0}^{n-1}a_j*(\sum_{i=0}^{n−1}(\omega_n^i)^{j-k}) j=0n1aj(i=0n1(ωni)jk)

j − k = 0 j−k=0 jk=0时, ∑ i = 0 n − 1 ( ω n i ) j − k = n \sum_{i=0}^{n−1}(\omega_n^i)^{j-k}=n i=0n1(ωni)jk=n
否则,通过等比数列求和可知:
∑ i = 0 n − 1 ( ω n i ) j − k \sum_{i=0}^{n−1}(\omega_n^i)^{j-k} i=0n1(ωni)jk= ( ω n j − k ) n − 1 ω n j − k − 1 \frac{(ω_n^{j−k})^n-1}{ω_n^{j−k}-1} ωnjk1(ωnjk)n1= ( ω n n ) j − k − 1 ω n j − k − 1 \frac{(ω_n^{n})^{j-k}-1}{ω_n^{j−k}-1} ωnjk1(ωnn)jk1= 1 − 1 ω n j − k − 1 \frac{1-1}{ω_n^{j−k}-1} ωnjk111=0
(因为 ω n n \omega_n^n ωnn= ω n 0 \omega_n^0 ωn0=1)
所以
z k z_k zk= n ∗ a k n*a_k nak
a k a_k ak= z k n \frac {z_k} n nzk ,得证。

怎么求单位根的倒数呢?
单位根的倒数其实就是它的共轭复数 。不明白的可以看看前面共轭复数的介绍
到现在你已经完全学会FFT了,但写递归还是可能会超时,所以我们需要优化

优化:迭代FFT

在进行FFT时,我们要把各个系数不断分组并放到两侧,一个系数原来的位置和最终的位置的规律如下。
初始位置: ω n 0 \omega_n^0 ωn0 ω n 1 \omega_n^1 ωn1 ω n 2 \omega_n^2 ωn2 ω n 3 \omega_n^3 ωn3 ω n 4 \omega_n^4 ωn4 ω n 5 \omega_n^5 ωn5 ω n 6 \omega_n^6 ωn6 ω n 7 \omega_n^7 ωn7
第一轮后: ω n 0 \omega_n^0 ωn0 ω n 2 \omega_n^2 ωn2 ω n 4 \omega_n^4 ωn4 ω n 6 \omega_n^6 ωn6| ω n 1 \omega_n^1 ωn1 ω n 3 \omega_n^3 ωn3 ω n 5 \omega_n^5 ωn5 ω n 7 \omega_n^7 ωn7
第二轮后: ω n 0 \omega_n^0 ωn0 ω n 4 \omega_n^4 ωn4| ω n 2 \omega_n^2 ωn2 ω n 6 \omega_n^6 ωn6| ω n 1 \omega_n^1 ωn1 ω n 5 \omega_n^5 ωn5| ω n 3 \omega_n^3 ωn3 ω n 7 \omega_n^7 ωn7
第三轮后: ω n 0 \omega_n^0 ωn0| ω n 4 \omega_n^4 ωn4| ω n 2 \omega_n^2 ωn2| ω n 6 \omega_n^6 ωn6| ω n 1 \omega_n^1 ωn1| ω n 5 \omega_n^5 ωn5| ω n 3 \omega_n^3 ωn3| ω n 7 \omega_n^7 ωn7
“|”代表分组界限
把每个位置用二进制表现出来。位置x上的数,最后所在的位置是“x二进制翻转得到的数”,例如4(100)最后到了1(001)5(101)最后不变为5(101),3(011)最后到了6(110)。
所以我们先把每个数放到最后的位置上,然后不断向上还原,同时求出点值表示就可以啦。
迭代版FFT就比之前的递归版快多了,真√ O ( n l o g 2 n ) O(nlog_2n) O(nlog2n)绝妙算法

代码实现FFT

下面是本人写的FFT加速高精度乘法的代码(并有详细注释):

#include<bits/stdc++.h>
using namespace std;
//complex是stl自带的定义复数的容器 
typedef complex<double> cp;
#define N 2097153
//pie表示圆周率π 
const double pie=acos(-1);
int n;
cp a[N],b[N];
int rev[N],ans[N];
char s1[N],s2[N];
//读入优化 
int read(){
    
	int sum=0,f=1;
	char ch=getchar();
	while(ch>'9'||ch<'0'){
    if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){
    sum=(sum<<3)+(sum<<1)+ch-'0';ch=getchar();}
	return sum*f;
}
//初始化每个位置最终到达的位置 
{
    
    int len=1<<k;
	for(int i=0;i<len;i++)
	rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}
//a表示要操作的系数,n表示序列长度
//若flag为1,则表示FFT,为-1则为IFFT(需要求倒数) 
void fft(cp *a,int n,int flag){
     
    for(int i=0;i<n;i++)
	{
    
	 //i小于rev[i]时才交换,防止同一个元素交换两次,回到它原来的位置。 
	  if(i<rev[i])swap(a[i],a[rev[i]]);
	}
	for(int h=1;h<n;h*=2)//h是准备合并序列的长度的二分之一
	{
    
	cp wn=exp(cp(0,flag*pie/h));//求单位根w_n^1 
	 for(int j=0;j<n;j+=h*2)//j表示合并到了哪一位
	 {
    
	  cp w(1,0);
	   for(int k=j;k<j+h;k++)//只扫左半部分,得到右半部分的答案
	   {
    
	     cp x=a[k];
	     cp y=w*a[k+h];
         a[k]=x+y;  //这两步是蝴蝶变换 
         a[k+h]=x-y;
         w*=wn; //求w_n^k 
	   }
	 }
	 }
	 //判断是否是FFT还是IFFT 
	 if(flag==-1)
	 for(int i=0;i<n;i++)
     a[i]/=n;
}
int main(){
    
	n=read(); 
	scanf("%s%s",s1,s2);
	//读入的数的每一位看成多项式的一项,保存在复数的实部 
    for(int i=0;i<n;i++)a[i]=(double)(s1[n-i-1]-'0');
	for(int i=0;i<n;i++)b[i]=(double)(s2[n-i-1]-'0');
	//k表示转化成二进制的位数 
	int k=1,s=2;
 	while((1<<k)<2*n-1)k++,s<<=1;
	init(k);
	//FFT 把a的系数表示转化为点值表示 
    fft(a,s,1);
    //FFT 把b的系数表示转化为点值表示 
    fft(b,s,1);
    //FFT 两个多项式的点值表示相乘 
    for(int i=0;i<s;i++)
    a[i]*=b[i];
    //IFFT 把这个点值表示转化为系数表示 
    fft(a,s,-1);
    //保存答案的每一位(注意进位) 
    for(int i=0;i<s;i++)
    {
    
    //取实数四舍五入,此时虚数部分应当为0或由于浮点误差接近0
	ans[i]+=(int)(a[i].real()+0.5);
	ans[i+1]+=ans[i]/10;
	ans[i]%=10;
	}
	while(!ans[s]&&s>-1)s--;
	if(s==-1)printf("0");
	else
	for(int i=s;i>=0;i--)
	printf("%d",ans[i]);
	return 0;
}

后记

这篇博客写了一天,终于写完了,完结撒花✿✿ヽ(°▽°)ノ✿
FWT我来啦!!!

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/Flag_z/article/details/99163939

智能推荐

获取大于等于一个整数的最小2次幂算法(HashMap#tableSizeFor)_整数 最小的2的几次方-程序员宅基地

文章浏览阅读2w次,点赞51次,收藏33次。一、需求给定一个整数,返回大于等于该整数的最小2次幂(2的乘方)。例: 输入 输出 -1 1 1 1 3 4 9 16 15 16二、分析当遇到这个需求的时候,我们可能会很容易想到一个"笨"办法:..._整数 最小的2的几次方

Linux 中 ss 命令的使用实例_ss@,,x,, 0-程序员宅基地

文章浏览阅读865次。选项,以防止命令将 IP 地址解析为主机名。如果只想在命令的输出中显示 unix套接字 连接,可以使用。不带任何选项,用来显示已建立连接的所有套接字的列表。如果只想在命令的输出中显示 tcp 连接,可以使用。如果只想在命令的输出中显示 udp 连接,可以使用。如果不想将ip地址解析为主机名称,可以使用。如果要取消命令输出中的标题行,可以使用。如果只想显示被侦听的套接字,可以使用。如果只想显示ipv4侦听的,可以使用。如果只想显示ipv6侦听的,可以使用。_ss@,,x,, 0

conda activate qiuqiu出现不存在activate_commandnotfounderror: 'activate-程序员宅基地

文章浏览阅读568次。CommandNotFoundError: 'activate'_commandnotfounderror: 'activate

Kafka 实战 - Windows10安装Kafka_win10安装部署kafka-程序员宅基地

文章浏览阅读426次,点赞10次,收藏19次。完成以上步骤后,您已在 Windows 10 上成功安装并验证了 Apache Kafka。在生产环境中,通常会将 Kafka 与外部 ZooKeeper 集群配合使用,并考虑配置安全、监控、持久化存储等高级特性。在生产者窗口中输入一些文本消息,然后按 Enter 发送。ZooKeeper 会在新窗口中运行。在另一个命令提示符窗口中,同样切换到 Kafka 的。Kafka 服务器将在新窗口中运行。在新的命令提示符窗口中,切换到 Kafka 的。,应显示已安装的 Java 版本信息。_win10安装部署kafka

【愚公系列】2023年12月 WEBGL专题-缓冲区对象_js 缓冲数据 new float32array-程序员宅基地

文章浏览阅读1.4w次。缓冲区对象(Buffer Object)是在OpenGL中用于存储和管理数据的一种机制。缓冲区对象可以存储各种类型的数据,例如顶点、纹理坐标、颜色等。在渲染过程中,缓冲区对象中存储的数据可以被复制到渲染管线的不同阶段中,例如顶点着色器、几何着色器和片段着色器等,以完成渲染操作。相比传统的CPU访问内存,缓冲区对象的数据存储和管理更加高效,能够提高OpenGL应用的性能表现。_js 缓冲数据 new float32array

四、数学建模之图与网络模型_图论与网络优化数学建模-程序员宅基地

文章浏览阅读912次。(1)图(Graph):图是数学和计算机科学中的一个抽象概念,它由一组节点(顶点)和连接这些节点的边组成。图可以是有向的(有方向的,边有箭头表示方向)或无向的(没有方向的,边没有箭头表示方向)。图用于表示各种关系,如社交网络、电路、地图、组织结构等。(2)网络(Network):网络是一个更广泛的概念,可以包括各种不同类型的连接元素,不仅仅是图中的节点和边。网络可以包括节点、边、连接线、路由器、服务器、通信协议等多种组成部分。网络的概念在各个领域都有应用,包括计算机网络、社交网络、电力网络、交通网络等。_图论与网络优化数学建模

随便推点

android 加载布局状态封装_adnroid加载数据转圈封装全屏转圈封装-程序员宅基地

文章浏览阅读1.5k次。我们经常会碰见 正在加载中,加载出错, “暂无商品”等一系列的相似的布局,因为我们有很多请求网络数据的页面,我们不可能每一个页面都写几个“正在加载中”等布局吧,这时候将这些状态的布局封装在一起就很有必要了。我们可以将这些封装为一个自定布局,然后每次操作该自定义类的方法就行了。 首先一般来说,从服务器拉去数据之前都是“正在加载”页面, 加载成功之后“正在加载”页面消失,展示数据;如果加载失败,就展示_adnroid加载数据转圈封装全屏转圈封装

阿里云服务器(Alibaba Cloud Linux 3)安装部署Mysql8-程序员宅基地

文章浏览阅读1.6k次,点赞23次,收藏29次。PS: 如果执行sudo grep 'temporary password' /var/log/mysqld.log 后没有报错,也没有任何结果显示,说明默认密码为空,可以直接进行下一步(后面设置密码时直接填写新密码就行)。3.(可选)当操作系统为Alibaba Cloud Linux 3时,执行如下命令,安装MySQL所需的库文件。下面示例中,将创建新的MySQL账号,用于远程访问MySQL。2.依次运行以下命令,创建远程登录MySQL的账号,并允许远程主机使用该账号访问MySQL。_alibaba cloud linux 3

excel离散度图表怎么算_excel离散数据表格-Excel 离散程度分析图表如何做-程序员宅基地

文章浏览阅读7.8k次。EXCEL中数据如何做离散性分析纠错。离散不是均值抄AVEDEV……=AVEDEV(A1:A100)算出来的是A1:A100的平均数。离散是指各项目间指标袭的离散均值(各数值的波动情况),数值较低表明项目间各指标波动幅百度小,数值高表明波动幅度较大。可以用excel中的离散公式为STDEV.P(即各指标平均离散)算出最终度离散度。excel表格函数求一组离散型数据,例如,几组C25的...用exc..._excel数据分析离散

学生时期学习资源同步-JavaSE理论知识-程序员宅基地

文章浏览阅读406次,点赞7次,收藏8次。i < 5){ //第3行。int count;System.out.println ("危险!System.out.println(”真”);System.out.println(”假”);System.out.print(“姓名:”);System.out.println("无匹配");System.out.println ("安全");

linux 性能测试磁盘状态监测:iostat监控学习,包含/proc/diskstats、/proc/stat简单了解-程序员宅基地

文章浏览阅读3.6k次。背景测试到性能、压力时,经常需要查看磁盘、网络、内存、cpu的性能值这里简单介绍下各个指标的含义一般磁盘比较关注的就是磁盘的iops,读写速度以及%util(看磁盘是否忙碌)CPU一般比较关注,idle 空闲,有时候也查看wait (如果wait特别大往往是io这边已经达到了瓶颈)iostatiostat uses the files below to create ..._/proc/diskstat

glReadPixels读取保存图片全黑_glreadpixels 全黑-程序员宅基地

文章浏览阅读2.4k次。问题:在Android上使用 glReadPixel 读取当前渲染数据,在若干机型(华为P9以及魅族某魅蓝手机)上读取数据失败,glGetError()没有抓到错误,但是获取到的数据有误,如果将获取到的数据保存成为图片,得到的图片为黑色。解决方法:glReadPixels实际上是从缓冲区中读取数据,如果使用了双缓冲区,则默认是从正在显示的缓冲(即前缓冲)中读取,而绘制工作是默认绘制到后缓..._glreadpixels 全黑