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$取模。

　　$n\leq 10^{18},k\leq 30000$

题解

　　毛爷爷论文题。

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

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

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

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

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

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

　　则不难列出转移方程：

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

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

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

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

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

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

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

$$dp_{x+y,i}=\sum_{j=0}^{i}2^{jy}\binom{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,0\cdots k}$。

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

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

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

　　听说有人用$FFT$+多项式求逆+多项式除法……等等算法，复杂度$O(k\log 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;
}