zoukankan      html  css  js  c++  java
  • FFT

    一、 多项式的表示法:

    1.系数表示法:

    一个(n-1)(n)项多项式(f(x))可以表示为(sumlimits_{i=1}^{n-1}a_i*x^i)

    于是可以用每一项的系数表示:(f(x)={a_0,a_1,...a_{n-1}})

    2.点值表示法:

    将多项式看成一个函数,那么找(n)个不同的数(x)代入,可以得到(n)个不同的(y)

    显然这(n)((x,y))可以唯一确定一个多项式,因为这相当于(n)个方程的线性方程组。

    写出来就是:(f(x)={(x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n-1},f(x_{n-1}))})

    二、多项式乘法:

    1.朴素算法:

    我们现在要求(h(x)=f(x)*g(x))

    我们分别考虑上两种表示法做乘法的复杂度:

    <1>系数表示法:我们要枚举(f(x))(g(x))的每一位系数相乘,复杂度显然为(O(n^2))
    <2>点值表示法:两个点值表示法的多项式相乘是(O(n))的!!!
    举个例子:
    (f(x)={(x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n-1},f(x_{n-1}))})
    (g(x)={(x_0,g(x_0)),(x_1,g(x_1)),...,(x_{n-1},g(x_{n-1}))})
    那么我们只需要(O(n))枚举每一位,便可得到:
    (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^2))的,从点值转系数还要(O(n^3))消元。

    系数表示法似乎无法优化,但是点值表示似乎有优化的空间。

    2.复数

    OI-WIKI,讲的很详细。

    3.单位复根

    1.定义:

    以下的(n)都是(2)的整次幂。

    我们在系数表示转点值表示时会找(n)个不同的(x_i)代入,对于每个(x_i),我们都要算出(x_i^0,x_i^1,...,x_i^{n-1}),这显然是(O(n^2))的。

    我们考虑取一些特殊的(x)优化这个过程,比如我们可以找(+1,-1,+i,-i),这些数的若干次方都是(1),但是只找这(4)个显然是不够的。

    现在我们在实轴和虚轴形成的坐标系中画一个单位圆,那么这个圆上的所有点都满足我们的要求,即若干次方后为(1)

    我们再将这个圆(n)等分,从((1,0))开始标号,分别标为(0,1,2,...,n-1)

    比如(n=8)

    (omega_n^k)表示(n)等分时标号为(k)的点代表的值,根据复数的运算(模长相乘,角相加)可以得到:(omega_n^k=(omega_n^1)^k)

    我们将(omega_n^1)称为(n)次单位根。

    根据三角函数的公式可以推出(omega_n^k)的公式:(omega_n^k=cosfrac{k}{n}2pi+isinfrac{k}{n}2pi)

    于是(omega_n^0,omega_n^1,..,omega_n^{n-1})就是我们要代入的(x_0,x_1,..,x_{n-1})

    2.性质:

    <1>(omega_n^k=omega_{2n}^{2k})

    证明只需代入公式。

    <2>(omega_n^{k+frac{n}{2}}=-omega_n^k)

    证明:
    (w_n^{frac{n}{2}}=cosfrac{n}{2}frac{2pi}{n}+isinfrac{n}{2}frac{2pi}{n})
    (=cospi+i*sinpi)
    (=-1)
    <3>(omega_n^n=omega_n^0=1)

    现在我们找到了一组特殊的(x),但是我们依旧要将每个(x_i)暴力代入算,所以复杂度仍然为(O(n^2))

    三、FFT

    上面那种方法叫做离散傅里叶变换((DFT)),我们考虑分治。

    假设我们有个多项式为(f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1})

    我们按照奇偶分类并将奇数类都提出一个(x)
    (f(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3...+a_{n-1}x^{n-1}))
    (=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2...+a_{n-1}x^{n-2}))

    此时我们定义另外两个多项式:
    (f1(x)=a_0+a_2x^1+a_4x^2+...+a_{n-2}x^{frac{n}{2}-1})
    (f2(x)=a_1+a_3x^1+...+a_{n-1}x^{frac{n}{2}-1})

    我们发现(f(x)=f1(x^2)+xf2(x^2))!!!!

    (k<frac{n}{2})

    (omega_n^k)代入:

    (f(omega_n^k)=f1((omega_n^k)^2)+omega_n^kf2((omega_n^k)^2))
    (=f1(omega_n^{2k})+omega_n^kf2(omega_n^{2k}))
    (=f1(omega_{frac{n}{2}}^k)+omega_n^kf2(omega_{frac{n}{2}}^{k}))

    (omega_n^{k+frac{n}{2}})代入:

    (f(omega_n^{k+frac{n}{2}})=f1(omega_n^{2k+n})+omega_n^{k+frac{n}{2}}f2(omega_n^{2k+n}))
    (=f1(omega_n^{2k}omega_n^n)-omega_n^kf2(omega_n^{2k}omega_n^n))
    (=f1(omega_n^{2k})+omega_n^kf2(omega_n^{2k}))
    (=f1(omega_{frac{n}{2}}^k)-omega_n^kf2(omega_{frac{n}{2}}^{k}))

    上面两个式子的结果除了符号不同以外完全相同,于是我们只要知道(f1(omega_{frac{n}{2}}^k),f2(omega_{frac{n}{2}}^{k})),我们就能求出
    (f(omega_n^k),f(omega_n^{k+frac{n}{2}}))

    于是就可以递归分治了,复杂度显然是(O(nlog_2n))的。

    四、IFFT

    我们用(FFT)求出了两个多项式乘积的点值表示((x_i,y_i)),我们还要将其转为系数表示。

    考虑怎么将点值表示快速转为系数表示:

    这里有个结论:
    一个多项式在分治的过程中乘上单位根的共轭复数,分治完的每一项除以n即为原多项式的每一项系数。

    我们设(c_k=sumlimits_{i=0}^{n-1}y_i(omega_n^{-k})^i)

    那么((c_i,w_n^{-i}))就是多项式(H(x)=y_0+y_1x+...+y_{n-1}x^{n-1})的点值表示。

    先推一波式子:
    (c_k=sumlimits_{i=0}^{n-1}y_i(omega_n^{-k})^i)
    (=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^i)^j)(omega_n^{-k})^i)
    (=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i)(omega_n^{-k})^i)
    (=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^j)^i(omega_n^{-k})^i))
    (=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^{j-k})^i))

    (S(x)=sumlimits_{i=0}^{n-1}x^i)

    (k!=0)时:
    代入(omega_n^k)可得:
    (S(omega_n^k)=1+omega_n^k+(omega_n^k)^2+...+(omega_n^k)^{n-1})
    两边同乘(omega_n^k)得:
    (omega_n^kS(omega_n^k)=omega_n^k+(omega_n^k)^2+...+(omega_n^k)^{n})
    两式相减可得:
    ((omega_n^k-1)S(omega_n^k)=(omega_n^k)^n-1)
    于是可得:
    (S(omega_n^k)=frac{(omega_n^k)^n-1}{omega_n^k-1})
    (=frac{(omega_n^n)^k-1}{omega_n^k-1})
    (=frac{1-1}{omega_n^k-1})
    (=0)
    (k=0)时:
    显然(S(omega_n^k)=n)

    回到之前的式子:
    (=sumlimits_{i=0}^{n-1}(sumlimits_{j=0}^{n-1}a_j(omega_n^{j-k})^i))

    (j!=k)时中间那个式子都是(0)
    于是:
    (c_k=na_k)
    (a_k=frac{c_k}{n})
    证完了。

    五、迭代FFT

    递归的常数太大,我们需要非递归的方法。

    先放张经典图:

    我们发现:
    每个位置分治后的最终位置为其二进制翻转后得到的位置。

    于是我们预先处理好每个数的最终位置,之后向上合并即可。

    模板题

    code:

    #include<bits/stdc++.h>
    using namespace std;
    const int maxn=1e7+10;
    const double Pi=acos(-1.0);
    int n,m,lim=1,len;
    int pos[maxn];
    struct cplx{double x,y;}a[maxn],b[maxn];
    cplx operator+(cplx a,cplx b){return (cplx){a.x+b.x,a.y+b.y};}
    cplx operator-(cplx a,cplx b){return (cplx){a.x-b.x,a.y-b.y};} 
    cplx operator*(cplx a,cplx b){return (cplx){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
    inline int read()
    {
    	char c=getchar();int res=0,f=1;
    	while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    	while(c>='0'&&c<='9')res=res*10+c-'0',c=getchar();
    	return res*f;
    }
    void fft(cplx* a,int op)
    {
    	for(int i=0;i<lim;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
    	for(int mid=1;mid<lim;mid<<=1)
    	{
    		cplx wn=(cplx){cos(Pi/mid),op*sin(Pi/mid)};
    		for(int i=0,r=mid<<1;i<lim;i+=r)
    		{
    			cplx w=(cplx){1,0};
    			for(int j=0;j<mid;j++,w=w*wn)
    			{
    				cplx x=a[i+j],y=w*a[i+mid+j];
    				a[i+j]=x+y;a[i+mid+j]=x-y;
    			}
    		}
    	}
    }
    int main()
    {
    	n=read(),m=read();
    	for(int i=0;i<=n;i++)a[i].x=read();
    	for(int i=0;i<=m;i++)b[i].x=read();
    	while(lim<=n+m)lim<<=1,len++;
    	for(int i=0;i<lim;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
    	fft(a,1);fft(b,1);
    	for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
    	fft(a,-1);
    	for(int i=0;i<=n+m;i++)printf("%d ",(int)(a[i].x/lim+0.5));
    	return 0;
    }
    
  • 相关阅读:
    java多线程详解(1)-多线程入门
    有关java中的hashCode问题
    java自动装箱和自动拆箱注意细节
    jquery选择器
    win10专业版激活密钥
    python3小例子:scrapy+mysql
    java List 等份截取
    七大查找算法
    十大经典排序算法
    jQuery.extend()方法
  • 原文地址:https://www.cnblogs.com/nofind/p/12118664.html
Copyright © 2011-2022 走看看