zoukankan      html  css  js  c++  java
  • 分治FFT学习笔记

    分治FFT m FFT与其优化


    前置知识 :

    快速傅里叶变换,快速数论变换,多项式求逆,CDQ分治,生成函数。
    【(蒟蒻的学习笔记)FFT m FFT NTT m NTT【或者可以自行去看其他更好的博客讲解】


    进入正题

    • 引入问题

    已知函数g(x)g(x)x[1,n1]xin [1,n-1]的函数值,求给定的f(x)f(x)的函数在x[1,n1]xin [1,n-1]的值,其中f(x)f(x)定义如下:
    f(x)=y=1xf(xy)g(y) f(x)=sumlimits_{y=1}^xf(x-y)g(y)
    其中,f(0)=1f(0)=1。(答案对998244353998244353取模)

    数据范围:2n1052leq nleq10^5,其余均小于模数。


    • 我会暴力!

    其实,根据定义式,直接就可以写出O(n2)O(n^2)(实际是n×(n1)2frac{n imes(n-1)}{2})的复杂度的算法。

    但是显然不行。

    • 我会FFT m FFT

    这个其实很容易看出是一个卷积式子,但是对于每一个f(x)f(x),我们不可能重新的去算一次FFT m FFT,因为那样的话复杂度会退化成O(n2logn)O(n^2logn),还不如暴力了,所以我们要考虑优化。

    如果我们已经知道了f(1)f(n2)f(1)sim f(frac{n}{2}),那么我们可以先计算出这些已知的对于未知部分f(n2+1)f(n)f(frac{n}{2}+1)sim f(n)的贡献,我们看原式,发现对于一个已知的f(x)f(x),对于后面的f(i)f(i)都只会有f(ij)g(i)[ij=x]f(i-j)g(i)[i-j=x]的贡献。

    所以对于一个当前的区间lr,mid=l+r2lsim r,mid=leftlfloorfrac{l+r}{2} ight floor我们已知f(l)f(mid)f(l)sim f(mid),那么它对于f(mid+1)f(r)f(mid+1)sim f(r)的贡献,可以用如下方法快速计算O((rl+1)log(rl+1))O((r-l+1)log(r-l+1)):

    我们令A(i)=f(i+l){i[0,midl]}A(i)=f(i+l){iin[0,mid-l]},再令B(i)=g(i){i[1,rl]}B(i)=g(i){iin[1,r-l]},那么我们再令C(i)=A(i)B(i)C(i)=A(i)igotimes B(i)(卷积,相当于C(x)=A(i)×B(j)[i+j=x]C(x)=sum A(i) imes B(j)[i+j=x]),那么对于f(mid+1)f(r)f(mid+1)sim f(r)的其中的f(x)f(x)贡献就为C(xl1)C(x-l-1),这个就是用FFT m FFT实现了。

    那么对于整个区间[1,n][1,n],我们分治下去,先计算[1,n+12][1,frac{n+1}{2}],这时你已经知道前半部分的值,然后FFT m FFT计算贡献,加到后半部分,然后去算[n+12,n][frac{n+1}{2},n]的值,这也就是CDQ m CDQ分治的过程。

    总的复杂度为O(nlog2n)O(nlog^2n)(总共递归lognlogn层,每层元素nn个,FFT m FFTO(nlogn)O(nlogn),所以总复杂度为O(分治×FFT)O( ext{分治} imes ext{FFT}))。

    Luogu 模板
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    #define ll long long
    using namespace std;
    const int M=4e5+10;
    const ll mod=998244353ll,G=3;//模数与原根
    ll inv[M],g[M],f[M];
    int R[M],lg,sze;
    
    ll fpow(ll a,ll b){
    	ll ans=1;
    	for(;b;b>>=1,a=(a*a)%mod){
    		if(b&1) ans=(ans*a)%mod;
    	}
    	return ans;
    }
    
    void NTT(ll *a,int n,int f){
    	for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);
    	for(int i=2;i<=n;i<<=1){
    		int now=i>>1;
    		ll wn=inv[i];//现场算或者预处理
    //		fpow(G,(mod-1)/i);
    		for(int j=0;j<n;j+=i){
    			ll w=1,x,y;
    			for(int k=j;k<j+now;k++,w=(w*wn)%mod){
    				x=a[k],y=w*a[k+now]%mod;
    				a[k]=(x+y)%mod;a[k+now]=((x-y)%mod+mod)%mod;
    			}
    		}
    	}
    	if(f==-1){
    		ll Inv=fpow(n,mod-2);
    		for(int i=0;i<=n;i++) a[i]=(a[i]*Inv)%mod;
    		for(int i=1;i<=(n>>1);i++)swap(a[i],a[n-i]);
    	}
    }
    
    void calc(ll *a,ll *b,int n){
    	NTT(a,n,1);NTT(b,n,1);
    	for(int i=0;i<=n;i++)a[i]=(a[i]*b[i])%mod;
    	NTT(a,n,-1);
    }
    
    ll A[M],B[M];
    void cdq(int l,int r){
    	if(l==r) return;
    	int mid=l+r>>1;
    	cdq(l,mid);
    
    	int up=r-l-1;
    //	for(lg=0,sze=1;sze<=(up<<1);sze<<=1)++lg;--lg;
    	for(lg=0,sze=1;sze<=up;sze<<=1)++lg;--lg;//其实只用的到区间长度,实际却有两倍,所以保险用两倍。
    	for(int i=0;i<sze;i++) R[i]=(R[i>>1]>>1)|((i&1)<<lg),A[i]=B[i]=0;
    //	memset(A,0,sizeof(A));memset(B,0,sizeof(B));
    //	fill(A,A+sze,0);fill(B,B+sze,0);//这两种方法赋值fill会比memset快一些(因为fill确定了大小)
    	for(int i=l;i<=mid;i++) A[i-l]=f[i];
    	for(int i=1;i<=r-l;i++) B[i-1]=g[i];
    	calc(A,B,sze);
    	
    	for(int i=mid+1;i<=r;i++) f[i]=(f[i]+A[i-l-1]%mod)%mod;
    	cdq(mid+1,r);
    }
    int n;
    int main(){
    	scanf("%d",&n);f[0]=1;
    	for(int i=1;i<n;i++)scanf("%lld",&g[i]);
    	for(int i=2,up=n<<2;i<=up;i<<=1)inv[i]=fpow(G,(mod-1)/i);
    	cdq(0,n-1);
    	for(int i=0;i<n;i++) printf("%lld%c",f[i],i==n-1?'
    ':' ');
    	return 0;
    }
    

    有没有更加优秀的做法呢?

    有,但是只有在取模的意义下,我们只需将式子变形即可。

    我们令FFff的生成函数,GGgg的生成函数。

    那么可以得知:

    F(x)=i=0f(i)xi F(x)=sumlimits_{i=0}^{infty}f(i)x^i

    G(x)=i=0g(i)xi G(x)=sumlimits_{i=0}^{infty}g(i)x^i

    那么求FFGG的卷积就有:

    F(x)G(x)=i=0xij+k=if(j)g(k) F(x)igotimes G(x)=sumlimits_{i=0}^{infty}x^isumlimits_{j+k=i}f(j)g(k)

    我们发现xij+k=if(j)g(k)x^isumlimits_{j+k=i}f(j)g(k)是不是和F(i)F(i)很像,所以我们可以知道那个卷积式可以写成F(x)f(0)x0F(x)-f(0)x^0,也就是F(x)f(0)F(x)-f(0),那么我们最后求的答案变为下面这个式子:
    F(x)G(x)F(x)f(0)(mod&ThinSpace;&ThinSpace;xn)F(x)f(0)1G(x)(mod&ThinSpace;&ThinSpace;xn)F(x)(1G(x))1(mod&ThinSpace;&ThinSpace;xn) F(x)igotimes G(x)equiv F(x)-f(0)(mod x^n) \ F(x)equiv frac{f(0)}{1-G(x)}(mod x^n) \ F(x)equiv (1-G(x))^{-1}(mod x^n)

    那么我们直接使用多项式求逆即可解决,复杂度O(nlogn)O(nlogn)


    博主也是才学习,有错误或者疑惑可以提出,大佬带带蒟蒻

  • 相关阅读:
    11,Django组件分页器
    10,Django于ajax
    阿里云安装Nexus搭建Maven私有仓库
    maven 自动部署到tomcat
    linux 7.2 下安装maven
    小程序防止遮罩层穿透
    Linux 下安装JDK
    Linux 命令未自动提示补全
    nginx 、tomcat 集群配置、shiro Session 共享
    nginx负载均衡配置
  • 原文地址:https://www.cnblogs.com/VictoryCzt/p/10053418.html
Copyright © 2011-2022 走看看