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;
    }
    
  • 相关阅读:
    tile38 复制配置
    The Guardian’s Migration from MongoDB to PostgreSQL on Amazon RDS
    tile38 一款开源的geo 数据库
    sqler sql 转rest api 的docker 镜像构建(续)使用源码编译
    sqler sql 转rest api javascript 试用
    sqler sql 转rest api redis 接口使用
    sqler sql 转rest api 的docker image
    sqler sql 转rest api 的工具试用
    apache geode 试用
    benthos v1 的一些新功能
  • 原文地址:https://www.cnblogs.com/nofind/p/12118664.html
Copyright © 2011-2022 走看看