zoukankan      html  css  js  c++  java
  • CodeForces 623E Transforming Sequence 动态规划 倍增 多项式 FFT 组合数学

    原文链接http://www.cnblogs.com/zhouzhendong/p/8848990.html

    题目传送门 - CodeForces 623E

    题意

      给定$n,k$。

      让你构造序列$a(0<a_i<2^k)$,满足$b_i(b_i=a_1 or a_2 or cdots or a_i)$严格单调递增。($or$为按位或)

      问你方案总数。对$10^9+7$取模。

      $nleq 10^{18},kleq 30000$

    题解

      毛爷爷论文题。

      我怀疑我看到的是假论文。里面不仅题面有点小问题,题解也有大问题。导致我一脸懵逼了好久QAQ。

      做这题时大概是觉得今天效率太低了,就$SKIP$了猜题解过程直接看题解了。

      考虑你有$k$个数位,要保证单调递增,由于是$or$运算,所以值为$1$的数位不管怎样或都不变了,所以每多一个$a_i$至少要多让一个数位变成$1$。

      因为只有$k$个数位,所以$a$的元素个数最多有$n$个,即$nleq k$。当$n>k$时就是无解,输出$0$。

      现在我们首先考虑大力$DP$。

      设$dp_{i,j}$表示前$i$个数占用了$j$个二进制位的方案数。

      则不难列出转移方程:

    $$dp_{i+1,x}=sum_{j=0}^{x-1}2^jinom{x}{j}dp_{i,j}$$

      注意$j$从$0$开始枚举其实和从$i$开始是等价的。

      因为原本就有$1$的数位,新增的数的那一位不管是$0$还是$1$都等价,所以要乘$2^j$。

      $inom{x}{j}$表示把$j$个二进制位$1$插入到新的$x$个二进制位$1$中的方案数。

      但是这样太慢了,要$TLE$。

      我们考虑升级一下上面的那个转移方程。

      设原本有$x$个数,现在一下子加入$y$个数,写出转移方程:

    $$dp_{x+y,i}=sum_{j=0}^{i}2^{jy}inom{i}{j}dp_{x,j}cdot dp_{y,i-j}$$

      其中,由于要加入$y$个数字,每一个数字都在原有的$j$个二进制位$1$中贡献了$2^j$的系数(和之前同理),所以总的系数贡献就是$(2^j)^y=2^{jy}$了。

      组合数意义还是和之前一样了。

      我们发现这个是个多项式卷积的形式,所以可以直接$FFT$优化。

      注意这个$FFT$要拆系数,不然要爆$LL$。

      这样的话,我们就不用一次一次转移,可以通过之前得到的结果来一坨一坨转移。

      然后我们发现,我们可以运用倍增快速幂的做法来解决整个运算过程。

      我们要求的是$dp_{n,0cdots k}$。

      最终的答案就是$sum_{i=0}^{k}inom{k}{i}dp_{n,i}$。

      注意,我代码里面把$k$写成了$m$。

      时间复杂度$O(klog^2 k)$。

      听说有人用$FFT$+多项式求逆+多项式除法……等等算法,复杂度$O(klog k)$干掉了此题??但是由于大常数却没快多少……

    代码

    #include <bits/stdc++.h>
    using namespace std;
    typedef long long LL;
    LL read(){
    	LL x=0;
    	char ch=getchar();
    	while (!('0'<=ch&&ch<='9'))
    		ch=getchar();
    	while ('0'<=ch&&ch<='9')
    		x=(x<<1)+(x<<3)+ch-48,ch=getchar();
    	return x; 
    }
    const int N=1<<17;
    const LL mod=1e9+7;
    double PI=acos(-1.0);
    int n,m;
    int s,d,R[N];
    LL Fac[N],Inv[N];
    LL Pow(LL x,LL y){
    	if (!y)
    		return 1LL;
    	LL xx=Pow(x,y/2);
    	xx=xx*xx%mod;
    	if (y&1LL)
    		xx=xx*x%mod;
    	return xx;
    }
    struct C{
    	double r,i;
    	C(){}
    	C(double a,double b){r=a,i=b;}
    	C operator + (C x){return C(r+x.r,i+x.i);}
    	C operator - (C x){return C(r-x.r,i-x.i);}
    	C operator * (C x){return C(r*x.r-i*x.i,r*x.i+i*x.r);}
    }w[N],A[N],B[N];
    LL tot[N],dp[N],now[N],D[N],E[N];
    void FFT(C a[],int n){
    	for (int i=0;i<n;i++)
    		if (i<R[i])
    			swap(a[i],a[R[i]]);
    	for (int t=n>>1,d=1;d<n;d<<=1,t>>=1)
    		for (int i=0;i<n;i+=(d<<1))
    			for (int j=0;j<d;j++){
    				C tmp=w[t*j]*a[i+j+d];
    				a[i+j+d]=a[i+j]-tmp;
    				a[i+j]=a[i+j]+tmp;
    			}
    }
    void FFT_mul(LL a[],LL b[],LL c[],int n){
    	for (int i=0;i<n;i++)
    		c[i]=0,A[i]=C(a[i]>>15,0),B[i]=C(b[i]>>15,0);
    	FFT(A,n),FFT(B,n);
    	for (int i=0;i<n;i++)
    		A[i]=A[i]*B[i],w[i].i*=-1.0;
    	FFT(A,n);
    	for (int i=0;i<n;i++){
    		tot[i]=((LL)(A[i].r/n+0.5))%mod;
    		c[i]=(c[i]+tot[i]*(1LL<<30))%mod;
    		w[i].i*=-1.0;
    	}
    	for (int i=0;i<n;i++)
    		A[i]=C(a[i]&32767,0),B[i]=C(b[i]&32767,0);
    	FFT(A,n),FFT(B,n);
    	for (int i=0;i<n;i++)
    		A[i]=A[i]*B[i],w[i].i*=-1.0;
    	FFT(A,n);
    	for (int i=0;i<n;i++){
    		tot[i]=(tot[i]+((LL)(A[i].r/n+0.5)))%mod;
    		c[i]=(c[i]+((LL)(A[i].r/n+0.5)))%mod;
    		w[i].i*=-1.0;
    	}
    	for (int i=0;i<n;i++){
    		A[i]=C((a[i]>>15)+(a[i]&32767),0);
    		B[i]=C((b[i]>>15)+(b[i]&32767),0);
    	}
    	FFT(A,n),FFT(B,n);
    	for (int i=0;i<n;i++)
    		A[i]=A[i]*B[i],w[i].i*=-1.0;
    	FFT(A,n);
    	for (int i=0;i<n;i++){
    		LL v=((LL)(A[i].r/n+0.5))%mod;
    		c[i]=(c[i]+(v-tot[i]+mod)%mod*(1LL<<15))%mod;
    		w[i].i*=-1.0;
    	}
    }
    void DP_mul(LL a[],LL b[],LL c[],int &x,int &y){
    	for (int i=0;i<s;i++)
    		D[i]=E[i]=0;
    	for (int i=0;i<=m;i++){
    		D[i]=a[i]*Inv[i]%mod*Pow(2,1LL*y*i%(mod-1))%mod;
    		E[i]=b[i]*Inv[i]%mod;
    	}
    	FFT_mul(D,E,c,s);
    	for (int i=0;i<s;i++)
    		c[i]=c[i]*(i<=m?Fac[i]:0)%mod;
    	x+=y;
    }
    void reads(){
    	LL nn=read();
    	scanf("%d",&m);
    	n=(nn>m)?-1:(int)nn;
    }
    int main(){
    	reads();
    	if (n==-1){
    		puts("0");
    		return 0;
    	}
    	Fac[0]=1;
    	for (int i=1;i<=m;i++)
    		Fac[i]=Fac[i-1]*i%mod;
    	Inv[m]=Pow(Fac[m],mod-2);
    	for (int i=m-1;i>=0;i--)
    		Inv[i]=Inv[i+1]*(i+1)%mod;
    	for (s=1,d=0;s<m*2+2;s<<=1,d++);
    	for (int i=0;i<s;i++){
    		R[i]=(R[i>>1]>>1)|((i&1)<<(d-1));
    		w[i]=C(cos(2*i*PI/s),sin(2*i*PI/s));
    	}
    	dp[0]=1;
    	for (int i=1;i<=m;i++)
    		now[i]=1;
    	int x,y;
    	for (x=0,y=1;y<=n;){
    		if (n&y)
    			DP_mul(dp,now,dp,x,y);
    		DP_mul(now,now,now,y,y);
    	}
    	LL ans=0;
    	for (int i=0;i<=m;i++)
    		ans=(ans+dp[i]*Fac[m]%mod*Inv[i]%mod*Inv[m-i])%mod;
    	printf("%I64d",ans);
    	return 0;
    }
    

      

  • 相关阅读:
    element ui 修改默认样式
    npm ERR! Cannot read property 'match' of undefined 错误处理
    Each record in table should have a unique `key` prop,or set `rowKey` to an unique primary key.
    (转)教你怎么理解正则表达式之零宽断言(环视)
    (转)通过扩展让ASP.NET Web API支持JSONP
    (转)走进AngularJs(六) 服务
    (转)在ASP.NET MVC3 中利用Jsonp跨域访问
    (转)DataTable与结构不同实体类之间的转换
    (转)深入理解最强桌面地图控件GMAP.NET --- 百度地图
    (转)打造一套UI与后台并重.net通用权限管理系统
  • 原文地址:https://www.cnblogs.com/zhouzhendong/p/CF623E.html
Copyright © 2011-2022 走看看