【CF623E】Transforming Sequence

by AmanoKumiko

Description

對於序列\(a\)，定義其變換\(b\)，其中\(b_j=\oplus_{i=1}^ja_i\)

給出\(n,k\)

計算有多少個長為\(n\)的序列\(a\)，滿足\(\forall a_i,1\le a_i\le 2^k-1\)，且\(b\)嚴格單調遞增

對\(1e9+7\)取模

Input

一行兩個數\(n,k\)

Output

一行一個數表示答案

Sample Input

3 3

Sample Output

48

Data Constraint

\(1\le n\le 10^{18},1\le k\le 3*10^4\)

Solution

首先由\(b\)的嚴格單增，可以發現每個\(a_i\)一定至少有一位滿足之前未出現，在位置\(i\)出現

那麼若\(n>k\)，答案就是\(0\)了

我們可以先選出要用的二進位制為\(1\)的位置的集合，再把集合分配給\(n\)個數

然後對於\(i\)，它分配的集合大小為\(S_i\)

設選出的集合大小為\(w\)，那麼總的貢獻就是\(\binom{k}{w}w!\prod_{i=1}^n\frac{2^{(n-i)|S_i|}}{|S_i|!}\)

這裡我們為了方便把\((n-i)\)變為\((i-1)\)，問題是等價的

於是變成了\(\sum_{i=n}^kk^{\underline i}2^{-i}\prod_{j=1}^n\frac{2^{j|S_j|}}{|S_j|!}\)

考慮設一個樸素的\(dp\)

那麼我們可以寫出生成函式

令\(F_i=\sum_{j=1}^{+\infty}\frac{2^{ij}x^j}{j!}\)

要求的就是\(\prod_{i=1}^nF_i\)

這裡給出一種思路不同於倍增FFT的做法

\[\begin{align} \prod_{i=1}^nF_i&=x^n2^{\frac{(1+n)n}{2}}\prod_{i=1}^n\frac{F_i}{2^ix}\\ &=x^n2^{\frac{(1+n)n}{2}}\exp \sum_{i=1}^nG(2^ix)\\ &=x^n2^{\frac{(1+n)n}{2}}\exp \sum_{j=0}^{+\infty}g_jx^j\sum_{i=1}^n2^{ij}\\ &=x^n2^{\frac{(1+n)n}{2}}\exp \sum_{j=0}^{+\infty}g_jx^j\frac{2^{(n+1)j}-2^j}{2^j-1}\\ \end{align} \]

\(G(x)=\ln \frac{e^{x}-1}{x}\)

那麼時間複雜度\(O(n\log n)\)

由於出題人比較毒瘤，給了\(1e9+7\)

所以還得寫個任意模數\(\exp\)。。。

Code

#include<bits/stdc++.h>
using namespace std;
#define F(i,a,b) for(int i=a;i<=b;i++)
#define Fd(i,a,b) for(int i=a;i>=b;i--)
#define mo 1000000007
#define N 600010
#define Ld long double

const Ld PI=acos(-1.0);

struct Complex{
	Ld real,img;
	Complex(Ld a=0.0,Ld b=0.0){real=a;img=b;}
	Complex Conj(){return Complex(real,-img);}
	friend Complex operator+(Complex a,Complex b){return Complex(a.real+b.real,a.img+b.img);}
	friend Complex operator-(Complex a,Complex b){return Complex(a.real-b.real,a.img-b.img);}
	friend Complex operator*(Complex a,Complex b){return Complex(a.real*b.real-a.img*b.img,a.real*b.img+a.img*b.real);}
	friend Complex operator*(Complex a,Ld b){return Complex(a.real*b,a.img*b);}
	friend Complex operator/(Complex a,Ld b){return Complex(a.real/b,a.img/b);}
	Complex operator+=(Complex x){return (*this)=(*this)+x;}
	Complex operator-=(Complex x){return (*this)=(*this)-x;}
	Complex operator*=(Complex x){return (*this)=(*this)*x;}
};

const Complex I=Complex(0,1);

int rev[N];

void BRT(int x){F(i,0,x-1)rev[i]=(rev[i>>1]>>1)|((i&1)?x>>1:0);}

struct poly_FFT{
	vector<Complex>val;
	void rsz(int x){val.resize(x);}
	int sz(){return val.size();}
	Complex operator[](int x)const{
		if(x<0||x>=val.size())return Complex();
		return val[x];
	}
	void FFT(int x){
		int n=sz();
		F(i,0,n-1)if(rev[i]<i)swap(val[rev[i]],val[i]);
		for(int mid=1;mid<n;mid<<=1){
			Complex wn=Complex(cos(PI/mid),sin(x*PI/mid));
			for(int i=0;i<n;i+=mid<<1){
				Complex w=Complex(1);
				F(j,0,mid-1){
					Complex t=w*val[i|j|mid];w*=wn;
					val[i|j|mid]=val[i|j]-t;val[i|j]+=t;
				}
			}
		}
		if(x==-1){F(i,0,n-1)val[i]=val[i]/(Ld)n;}
	}
	void DFT(){FFT(1);}
	void IDFT(){FFT(-1);}
};


int num(Ld x){
	long long y=x<0?x-0.5:x+0.5;
	return (y%mo+mo)%mo;
}

const int M=sqrt(mo)+1;

void MTT(int*A,int*B,int*C,int l){
	poly_FFT Pa,Pb,a0,a1,b0,b1;
	Pa.rsz(l);Pb.rsz(l);a0.rsz(l);a1.rsz(l);b0.rsz(l);b1.rsz(l);
	F(i,0,l-1)Pa.val[i]=Complex(A[i]/M,A[i]%M),Pb.val[i]=Complex(B[i]/M,B[i]%M);
	Pa.DFT();Pb.DFT();
	F(i,0,l-1){
		Complex Qa=(i?Pa[l-i].Conj():Pa[0].Conj()),Qb=(i?Pb[l-i].Conj():Pb[0].Conj());
		a0.val[i]=(Pa[i]+Qa)*0.5;a1.val[i]=(Qa-Pa[i])*0.5*I;b0.val[i]=(Pb[i]+Qb)*0.5;b1.val[i]=(Qb-Pb[i])*0.5*I;
	}
	F(i,0,l-1)Pa.val[i]=a0[i]*b0[i]+I*a1[i]*b0[i],Pb.val[i]=a0[i]*b1[i]+I*a1[i]*b1[i];
	Pa.IDFT();Pb.IDFT();
	F(i,0,l-1)C[i]=(1ll*M*M%mo*num(Pa[i].real)%mo+1ll*M*(num(Pa[i].img)+num(Pb[i].real))%mo+1ll*num(Pb[i].img)%mo)%mo;
}

int G1[N],G2[N],fac[N],ifac[N],inv[N];

int mod(int x){return x>=mo?x-mo:x;}

int mi(int x,int y){
	if(!y)return 1;
	if(y==1)return x;
	return y%2?1ll*x*mi(1ll*x*x%mo,y/2)%mo:mi(1ll*x*x%mo,y/2);
}

void init(){
	fac[0]=ifac[0]=1;
	F(i,1,N-10)fac[i]=1ll*fac[i-1]*i%mo,inv[i]=(i==1?1:1ll*mo/i*mod(mo-1ll*inv[mo%i]%mo)%mo);
	ifac[N-10]=mi(fac[N-10],mo-2);
	Fd(i,N-11,1)ifac[i]=1ll*ifac[i+1]*(i+1)%mo;
}

struct poly_NTT{
	vector<int>val;
	poly_NTT(int x=0){if(x)val.push_back(x);}
	poly_NTT(const vector<int>&x){val=x;}
	void Rev(){reverse(val.begin(),val.end());}
	void ins(int x){val.push_back(x);}
	void clear(){vector<int>().swap(val);}
	int sz(){return val.size();}
	void rsz(int x){val.resize(x);}
	void shrink(){for(;sz()&&!val.back();val.pop_back());}
	poly_NTT modxn(int x){
		if(val.size()<=x)return poly_NTT(val);
		else return poly_NTT(vector<int>(val.begin(),val.begin()+x));
	}
	int operator[](int x)const{
		if(x<0||x>=val.size())return 0;
		return val[x];
	}
	friend poly_NTT operator*(poly_NTT x,poly_NTT y){
		if(x.sz()<30||y.sz()<30){
			if(x.sz()>y.sz())swap(x,y);
			poly_NTT ret;
			ret.rsz(x.sz()+y.sz());
			F(i,0,ret.sz()-1){
				for(int j=0;j<=i&&j<x.sz();j++)
					ret.val[i]=mod(ret.val[i]+1ll*x[j]*y[i-j]%mo);
			}
	//		ret.shrink();
			return ret;
		}
		int l=1;
		while(l<x.sz()+y.sz()-1)l<<=1;
		x.rsz(l);BRT(l);
		static int A[N],B[N],C[N];
		F(i,0,x.sz()-1)A[i]=x[i];F(i,x.sz(),l-1)A[i]=0;
		F(i,0,y.sz()-1)B[i]=y[i];F(i,y.sz(),l-1)B[i]=0;
		MTT(A,B,C,l);
		F(i,0,l-1)x.val[i]=C[i];
//		x.shrink();
		return x;
	}
	friend poly_NTT operator+(poly_NTT x,poly_NTT y){
		poly_NTT ret;
		ret.rsz(max(x.sz(),y.sz()));
		F(i,0,ret.sz()-1)ret.val[i]=mod(x[i]+y[i]);
		return ret;
	}
	friend poly_NTT operator-(poly_NTT x,poly_NTT y){
		poly_NTT ret;
		ret.rsz(max(x.sz(),y.sz()));
		F(i,0,ret.sz()-1)ret.val[i]=mod(x[i]-y[i]+mo);
		return ret;
	}
	poly_NTT &operator*=(poly_NTT x){return (*this)=(*this)*x;}
	poly_NTT &operator+=(poly_NTT x){return (*this)=(*this)+x;}
	poly_NTT &operator-=(poly_NTT x){return (*this)=(*this)-x;}
	poly_NTT deriv(){
		poly_NTT f;
		f.rsz(sz()-1);
		F(i,0,sz()-2)f.val[i]=1ll*(i+1)*val[i+1]%mo;
		return f;
	}
	poly_NTT integ(){
		poly_NTT f;
		f.rsz(sz()+1);
		F(i,1,sz())f.val[i]=1ll*val[i-1]*inv[i]%mo;
		return f;
	}
	poly_NTT inver(int Len){
		poly_NTT res(mi(val[0],mo-2));
		for(int i=1;i<Len;){
			i<<=1;
			res=(res*(2-res*modxn(i))).modxn(i);
		}
		return res.modxn(Len);
	}
	poly_NTT Sqrt(int Len){
		poly_NTT f,res(1);
		for(int i=1;i<Len;){
			i<<=1;
			f=res.inver(i);
			res=((res+f*modxn(i))*inv[2]).modxn(i);
		}
		return res.modxn(Len);
	}
	poly_NTT Ln(int Len){
		return (deriv()*inver(Len)).integ().modxn(Len);
	}
	poly_NTT Exp(int Len){
		poly_NTT f(1);
		for(int i=2;i<Len*2;i<<=1)f=(f*(1-f.Ln(i)+modxn(i))).modxn(i);
		return f.modxn(Len);
	}
};

int n,k,ans;

int main(){
	init();
	long long x; 
	scanf("%lld",&x);
	scanf("%d",&k);
	if(x>k){printf("0");return 0;}
	n=x;
	poly_NTT g;
	F(i,0,k-n)g.ins(ifac[i+1]);
	g=g.Ln(k-n+1);
	int tmp=1;
	F(i,0,k-n){
		int v=1ll*mod(mi(2,(n+1)*i)-tmp+mo)*mi(tmp-1,mo-2)%mo;
		tmp=1ll*tmp*2%mo;
		g.val[i]=1ll*g[i]*v%mo;
	}
	g=g.Exp(k-n+1);
	tmp=mi(2,(1+n)*n/2);
	F(i,0,k-n)g.val[i]=1ll*g[i]*tmp%mo;
	tmp=mi(mi(2,n),mo-2);
	F(i,n,k){
		ans=mod(ans+1ll*g[i-n]*fac[k]%mo*ifac[k-i]%mo*tmp%mo),
		tmp=1ll*tmp*inv[2]%mo;
	}
	printf("%d",ans);
	return 0;
}