精华内容
下载资源
问答
  • 十分简明易懂的FFT(快速傅里叶变换

    万次阅读 多人点赞 2018-08-07 11:38:56
    快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶...

    FFT前言

    快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。

    FFT(Fast Fourier Transformation) 是离散傅氏变换(DFT)的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。

    ——百度百科

    FFT(Fast Fourier Transformation),中文名快速傅里叶变换,用来加速多项式乘法
    朴素高精度乘法时间O(n2)O(n^2),但FFTFFTO(nlog2n)O(n\log_2 n)的时间解决
    FFTFFT名字逼格高,也难懂,其他教程写得让人看不太懂,于是自己随便写一下

    • 建议对复数三角函数相关知识有所耳闻 (不会也无所谓)

    • 下面难懂的点我会从网上盗


    多项式的系数表示法和点值表示法

    • FFTFFT其实是一个用O(nlog2n)O(n\log_2n)的时间将一个用系数表示的多项式转换成它的点值表示的算法

    • 多项式的系数表示和点值表示可以互相转换

    系数表示法

    一个n-1次n项多项式f(x)f(x)可以表示为f(x)=i=0n1aixif(x)=\sum^{n-1}_{i=0}a_ix^i
    也可以用每一项的系数来表示f(x)f(x),即f(x)={a0,a1,a2,...,an1}f(x)=\{a_0,a_1,a_2,...,a_{n-1} \}
    这就是系数表示法,也就是平时数学课上用的方法

    点值表示法

    • 把多项式放到平面直角坐标系里面,看成一个函数

    • nn个不同的xx代入,会得出nn个不同的yy,在坐标系内就是nn个不同的点

    • 那么这nn个点唯一确定该多项式,也就是有且仅有一个多项式满足k,f(xk)=yk∀k,f(x_k)=y_k

    • 理由很简单,把nn条式子联立起来成为一个有n条方程的n元方程组,每一项的系数都可以解出来

    那么f(x)f(x)还可以用f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),...,(xn1,f(xn1))}f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\}来表示
    这就是点值表示法


    高精度乘法下两种多项式表示法的区别

    对于两个用系数表示的多项式,我们把它们相乘
    设两个多项式分别为A(x),B(x)A(x),B(x)
    我们要枚举AA每一位的系数与BB每一位的系数相乘
    那么系数表示法做多项式乘法时间复杂度O(n2)O(n^2)

    但两个用点值表示的多项式相乘,只需要O(n)O(n)的时间

    什么意思呢?

    设两个点值多项式分别为f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),...,(xn1,f(xn1))}f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\}g(x)={(x0,g(x0)),(x1,g(x1)),(x2,g(x2)),...,(xn1,g(xn1))}g(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),...,(x_{n-1},g(x_{n-1}))\}
    设它们的乘积是h(x)h(x),那么h(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),...,(xn1,f(xn1)g(xn1))}h(x)=\{(x_0,f(x_0)·g(x_0)),(x_1,f(x_1)·g(x_1)),...,(x_{n-1},f(x_{n-1})·g(x_{n-1}))\}

    所以这里的时间复杂度只有一个枚举的O(n)O(n)

    • 突然感觉高精度乘法能O(n)O(n)暴艹一堆题?

    • 但是朴素的系数表示法转点值表示法的算法还是O(n2)O(n^2)的,逆操作类似

    • 朴素系数转点值的算法叫DFT(离散傅里叶变换),点值转系数叫IDFT(离散傅里叶逆变换)

    • 难道高精度乘法只能O(n2)O(n^2)了吗?


    DFT前置知识&技能

    复数

    毕竟高中有所以不多说

    我们把形如a+bi(a,b均为实数)的数称为复数,其中a称为实部,b称为虚部,i称为虚数单位。当虚部等于零时,这个复数可以视为实数;当z的虚部不等于零时,实部等于零时,常称z为纯虚数。复数域是实数域的代数闭包,也即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入,经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作,此概念逐渐为数学家所接受。
    ——百度百科

    初中数学老师会告诉你没有1\sqrt{-1},但仅限RR
    扩展至复数集CC,定义i2=1i^2=-1,一个复数zz可以表示为z=a+bi(a,bR)z=a+bi(a,b\in R)
    其中aa称为实部bb称为虚部ii称为虚数单位

    • 在复数集中就可以用ii表示负数的平方根,如7=7i\sqrt{-7}=\sqrt{7}i

    还可以把复数看成平面直角坐标系上的一个点,比如下面


    xx轴就是实数集中的坐标轴,yy轴就是虚数单位ii

    这个点(2,3)(2,3)表示的复数就是2+3i2+3i,或者想象它代表的向量(2,3)(2,3)
    其实我们还可以自己想象 (其实没有这种表达方式) 它可以表示为(13,θ)(\sqrt{13},\theta)
    一个复数zz定义为它到原点的距离,记为z=a2+b2|z|=\sqrt{a^2+b^2}
    一个复数z=a+biz=a+bi共轭复数abia-bi(虚部取反),记为z=abi\overline{z}=a-bi

    复数的运算

    复数不像点或向量,它和实数一样可以进行四则运算
    设两个复数分别为z1=a+bi,z2=c+diz_1=a+bi,z_2=c+di,那么
    z1+z2=(a+c)+(b+d)iz_1+z_2=(a+c)+(b+d)iz1z2=(acbd)+(ad+bc)iz_1z_2=(ac−bd)+(ad+bc)i

    复数相加也满足平行四边形法则


    这张是从网上盗的

    AB+AD=ACAB+AD=AC

    复数相乘还有一个值得注意的小性质
    (a1,θ1)(a2,θ2)=(a1a2,θ1+θ2)(a_1,\theta_1)*(a_2,\theta_2)=(a_1a_2,\theta_1+\theta_2)
    模长相乘,极角相加


    DFT(离散傅里叶变换)

    • 一定注意从这里开始所有的nn都默认为22的整数次幂

    对于任意系数多项式转点值,当然可以随便取任意nnxx值代入计算
    但是暴力计算xk0,xk1,...,xkn1(k[0,n))x_k^0,x_k^1,...,x_k^{n-1}(k\in[0,n))当然是O(n2)O(n^2)的时间
    其实可以代入一组神奇xx,代入以后不用做那么多的次方运算
    这些xx当然不是乱取的,而且取这些xx值应该就是 傅里叶 的主意了

    考虑一下,如果我们代入一些xx,使每个xx若干次方等于11,我们就不用做全部的次方运算了
    ±1±1是可以的,考虑虚数的话±i±i也可以,但只有这四个数远远不够

    • 傅里叶说:这个圆圈上面的点都可以做到

    以原点为圆心,画一个半径为11单位圆
    那么单位圆上所有的点都可以经过若干次次方得到11
    傅里叶说还要把它给nn等分了,比如此时n=8n=8

    橙色点即为n=8n=8时要取的点,从(1,0)(1,0)点开始,逆时针从00号开始标号,标到77
    记编号为kk的点代表的复数值为ωnk\omega_n^k,那么由模长相乘,极角相加可知(ωn1)k=ωnk(\omega_n^1)^k=\omega_n^k
    其中ωn1\omega_n^1称为nn次单位根,而且每一个ω\omega都可以求出 (我三角函数不好)
    ωnk=coskn2π+isinkn2π\omega_n^k=\cos{k\over n}2π+i\sin{k\over n} 2π

    或者说也可以这样解释这条式子


    注意sin2θ+cos2θ=1sin^2\theta+cos^2\theta=1什么的,就容易理解了

    那么ωn0,ωn1,...,ωnn1\omega^0_n,\omega^1_n,...,\omega^{n-1}_n即为我们要代入的x0,x1,...,xn1x_0,x_1,...,x_{n-1}


    单位根的一些性质

    FFTFFT的过程中需要用到ω\omega的一些性质

    ωnk=ω2n2k\omega^k_n=\omega^{2k}_{2n}

    • 它们表示的点(或向量)表示的复数是相同的

    • 证明

    • ωnk=coskn2π+isinkn2π=cos2k2n2π+isin2k2n2π=ω2n2k\omega^k_n=cos{k\over n}2π+isin{k\over n} 2π=cos{2k\over 2n}2π+isin{2k\over 2n} 2π=\omega^{2k}_{2n}

    ωnk+n2=ωnk\omega^{k+{n \over 2}}_n=-\omega_n^k

    • 它们表示的点关于原点对称,所表示的复数实部相反,所表示的向量等大反向

    • 证明

    • ωnn2=cosn2n2π+isinn2n2π=cosπ+isinπ=1\omega^{n\over 2}_n=cos{{n\over 2}\over n}2\pi+isin{{n\over 2}\over n}2\pi=cos\pi+isin\pi=-1

    • (这个东西和eix=cosx+isinxe^{ix}=cosx+isinxeiπ+1=0e^{i\pi}+1=0有点关系,我不会就不讲了

    ωn0=ωnn\omega^0_n=\omega^n_n

    • 都等于11,或1+0i1+0i

    FFT(快速傅里叶变换)

    虽然DFTDFT搞出来一堆很牛逼ω\omega作为代入多项式的xx
    但只是代入这类特殊xx值法的变换叫做DFTDFT而已,还是要代入单位根暴力计算

    • DFT还是暴力O(n2)O(n^2)

    DFTDFT可以分治来做,于是 FFT(快速傅里叶变换) 就出来了
    首先设一个多项式A(x)A(x)
    A(x)=i=0n1aixi=a0+a1x+a2x2+...+an1xn1A(x)=\sum^{n-1}_{i=0}a_ix^i=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}

    A(x)A(x)下标的奇偶性A(x)A(x)分成两半,右边再提一个xx

    A(x)=(a0+a2x2+...+an2xn2)+(a1x+a3x3+...+an1xn1)A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})

    =(a0+a2x2+...+an2xn2)+x(a1+a3x2+...+an1xn2)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})

    两边好像非常相似,于是再设两个多项式A1(x),A2(x)A_1(x),A_2(x),令

    A1(x)=a0+a2x+a4x2+...+an2xn21A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{{n\over 2}-1}A2(x)=a1+a3x+a5x2+...+an1xn21A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{{n \over 2}-1}

    比较明显得出
    A(x)=A1(x2)+xA2(x2)A(x)=A_1(x^2)+xA_2(x^2)

    再设k<n2k<{n\over 2},把ωnk\omega^k_n作为xx代入A(x)A(x)(接下来几步变换要多想想单位根的性质)

    A(ωnk)=A1((ωnk)2)+ωnkA2((ωnk)2)A(\omega^k_n)=A_1((\omega^k_n)^2)+\omega^k_nA_2((\omega^k_n)^2)=A1(ωn2k)+ωnkA2(ωn2k)=A1(ωn2k)+ωnkA2(ωn2k)=A_1(\omega^{2k}_n)+\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})+\omega^k_nA_2(\omega^k_{n\over 2})

    那么对于A(ωnk+n2)A(\omega^{k+{n\over2}}_n)的话,代入ωnk+n2\omega^{k+{n \over 2}}_n
    A(ωnk+n2)=A1(ωn2k+n)+ωnk+n2A2(ωn2k+n)A(\omega^{k+{n\over 2}}_n)=A_1(\omega^{2k+n}_n)+\omega^{k+{n\over 2}}_nA_2(\omega^{2k+n}_n)=A1(ωn2kωnn)ωnkA2(ωn2kωnn)=A_1(\omega^{2k}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^n_n)=A1(ωn2k)ωnkA2(ωn2k)=A1(ωn2k)ωnkA2(ωn2k)=A_1(\omega^{2k}_n)-\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})-\omega^k_nA_2(\omega^k_{n\over2})

    • 发现了什么?

    A(ωnk)A(\omega^k_n)A(ωnk+n2)A(\omega^{k+{n\over2}}_n)两个多项式后面一坨东西只有符号不同
    就是说,如果已知A1(ωn2k)A_1(\omega^k_{n\over 2})A2(ωn2k)A_2(\omega^k_{n\over 2})的值,我们就可以同时知道A(ωnk)A(\omega^k_n)A(ωnk+n2)A(\omega^{k+{n\over2}}_n)的值
    现在我们就可以递归分治来搞FFTFFT

    每一次回溯时只扫当前前面一半的序列,即可得出后面一半序列的答案
    n==1n==1时只有一个常数项,直接returnreturn
    时间复杂度O(nlog2n)O(n\log_2n)


    IFFT(快速傅里叶逆变换)

    想一下,我们不仅要会FFTFFT,还要会IFFT(快速傅里叶逆变换)
    我们把两个多项式相乘 (也叫求卷积),做完两遍FFTFFT也知道了积的多项式的点值表示
    可我们平时用系数表示的多项式,点值表示没有意义

    • 怎么把点值表示的多项式快速转回系数表示法?

    • IDFTIDFT暴力O(n2)O(n^2)做?其实也可以用FFTFFTO(nlog2n)O(n\log_2n)的时间搞

    你有没有想过为什么傅里叶是把ωnk\omega^k_n作为xx代入而不是别的什么数?
    原因是有的但是有我也看不懂
    由于我是沙雕所以只用记住一个结论

    • 一个多项式在分治的过程中乘上单位根的共轭复数,分治完的每一项除以nn即为原多项式的每一项系数

    意思就是说FFTFFTIFFTIFFT可以一起搞


    朴素版FFT板子

    c++c++有自带的复数模板complexcomplex
    a.real()a.real()即表示复数aa的实部

    #include<complex>
    #define cp complex<double>
    
    void fft(cp *a,int n,int inv)//inv是取共轭复数的符号
    {
        if (n==1)return;
        int mid=n/2;
        static cp b[MAXN];
        fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1];
        fo(i,0,n-1)a[i]=b[i];
        fft(a,mid,inv),fft(a+mid,mid,inv);//分治
        fo(i,0,mid-1)
        {
            cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));//inv取决是否取共轭复数
            b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid];
        }
        fo(i,0,n-1)a[i]=b[i];
    }
    

    两个多项式a,ba,b相乘再转系数多项式cc,通常只用打这么一小段

    	cp a[MAXN],b[MAXN];
    	int c[MAXN];
    	fft(a,n,1),fft(b,n,1);//1系数转点值
    	fo(i,0,n-1)a[i]*=b[i];
    	fft(a,n,-1);//-1点值转系数
    	fo(i,0,n-1)c[i]=(int)(a[i].real()/n+0.5);//注意精度
    

    很明显,FFT只能处理nn22的整数次幂的多项式
    所以在FFTFFT前一定要把nn调到22的次幂

    这个板子看着好像很优美,但是

    递归常数太大,要考虑优化…


    FFTの优化——迭代版FFT

    这个图也是盗的

    这个很容易发现点什么吧?

    • 每个位置分治后的最终位置为其二进制翻转后得到的位置

    这样的话我们可以先把原序列变换好,把每个数放在最终的位置上,再一步一步向上合并
    一句话就可以O(n)O(n)预处理第ii位最终的位置rev[i]rev[i]

    fo(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    

    至于蝴蝶变换它死了其实是我不会


    真·FFT板子

    void fft(cp *a,int n,int inv)
    {
        int bit=0;
        while ((1<<bit)<n)bit++;
        fo(i,0,n-1)
        {
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
            if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次(就是没交换)
        }
        for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一
        {
        	cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根,pi的系数2已经约掉了
            for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度,i是合并到了哪一位
    		{
                cp omega(1,0);
                for (int j=0;j<mid;j++,omega*=temp)//只扫左半部分,得到右半部分的答案
                {
                    cp x=a[i+j],y=omega*a[i+j+mid];
                    a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换什么的
                }
            }
        }
    }
    

    这个板子好像不是那么好背
    至少这个板子已经很优美


    FFT后记

    本人版权意识薄弱……

    NTTNTT我来了

    展开全文
  • 快速傅里叶变换

    2019-01-20 22:31:18
    这是快速傅里叶变换的资料,可以通过该资料学习快速傅里叶变换
  • 精通matlab傅里叶变换及逆变换和快速傅里叶变换
  • 超详细易懂FFT(快速傅里叶变换)及代码实现

    万次阅读 多人点赞 2019-08-11 18:29:24
    前言 昨天学了一晚上,终于搞懂了FFT。希望能写一篇清楚易懂的题解分享给大家,也进一步加深自己的...FFT(Fast Fourier Transformation),中文名快速傅里叶变换,是离散傅氏变换的快速算法,它是根据离散傅氏变...

    前言

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

    定义

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

    有关知识

    多项式的两种表示方法

    系数表示法

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

    点值表示法

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

    多项式相乘

    设两个多项式分别为f(x)f(x),g(x)g(x),我们要把这两个多项式相乘 (即求卷积)。
    如果用系数表示法:
    我们要枚举ff的每一位的系数与gg的每一位的系数相乘,多项式乘法时间复杂度O(n2)O(n^2),这也是我们所熟知的高精度乘法的原理。
    如果用点值表示法:
    f[x]f[x]={(x0,f[x0])(x_0,f[x_0]),(x1,f[x1])(x_1,f[x_1]),…,(xn,f[xn])(x_n,f[x_n])}
    g[x]g[x]={(x0,g[x0])(x_0,g[x_0]),(x1,g[x1])(x_1,g[x_1]),…,(xn,g[xn])(x_n,g[x_n])}
    f[x]g[x]f[x]*g[x]={(x0,f[x0]g[x0])(x_0,f[x_0]*g[x_0]),(x1,f[x1]g[x1])(x_1,f[x_1]*g[x_1]),…,(xn,f[xn]g[xn])(x_n,f[x_n]*g[x_n])}
    我们可以发现,如果两个多项式取相同的xx,得到不同的yy值,那么只需要yy值对应相乘就可以了!
    复杂度只有枚举xxO(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表示的是傅里叶变换

    复数

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

    1.定义

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

    2.复数的几何意义

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

    3.运算法则

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

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

    4. 共轭复数

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

    FFT加速多项式乘法

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

    单位根

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

    怎么找单位根?

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

    单位根的性质

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

    1.ωnk=ω2n2k\omega_n^k=\omega_{2n}^{2k}

    2.ωnk=ωnk+n2\omega_n^k=-\omega_{n}^{k+\frac n 2}

    3.ωn0=ωnn=1\omega_n^0=\omega_{n}^n=1

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

    DFT(离散傅里叶变换)

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

    FFT(快速傅里叶变换)

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

    A(ωnk+n2)=A0(ωn2k+n)+ωnk+n2A1(ωn2k+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})
    =A0(ωn2k)ωnkA1(ωn2k)=A_0(\omega_\frac n2^{k})-\omega_n^{k}*A_1(\omega_\frac n 2^{k})

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

    IFFT(快速傅里叶逆变换)

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

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

    jk=0j−k=0时,i=0n1(ωni)jk=n\sum_{i=0}^{n−1}(\omega_n^i)^{j-k}=n
    否则,通过等比数列求和可知:
    i=0n1(ωni)jk\sum_{i=0}^{n−1}(\omega_n^i)^{j-k}=(ωnjk)n1ωnjk1\frac{(ω_n^{j−k})^n-1}{ω_n^{j−k}-1}=(ωnn)jk1ωnjk1\frac{(ω_n^{n})^{j-k}-1}{ω_n^{j−k}-1}=11ωnjk1\frac{1-1}{ω_n^{j−k}-1}=0
    (因为ωnn\omega_n^n=ωn0\omega_n^0=1)
    所以
    zkz_k=nakn*a_k
    aka_k=zkn\frac {z_k} n ,得证。

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

    优化:迭代FFT

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

    代码实现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我来啦!!!

    展开全文

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,493
精华内容 1,797
关键字:

快速傅里叶变换