首页 > 其他分享> 文章详细

【LuoguP5850】calc加强版

2022-01-06 19:00:09 阅读：177 来源： 互联网

标签：return 加强版 val int mo LuoguP5850 1ll calc sum

【LuoguP5850】calc加强版

by AmanoKumiko

Description

对于\(n=1,2,3 ... m\)，求出所有合法序列的权值

其中，合法的定义为每个元素都为\([1,k]\)中的整数，同时每个元素都不相同

同时，权值为所有元素的乘积

两个序列不同当且仅当他们任意一位不一样。

对\(998244353\)取模

Input

一行两个数\(k,m\)

Output

共\(m\)行表示答案

Sample Input

13 8

Sample Output

Data Constraint

\(1\le m\le k\le998244353\)

Solution

不难发现对于一个\(n\)，要求的就是\(n![x^n]\prod_{i=1}^k(ix+1)\)

由于\(k\)很大，所以不能分治\(NTT\)

但是可以运用积化和的\(trick\)

令\(F(x)=\prod_{i=1}^k(ix+1)\)

那么\(ln\ F(x)=\sum_{i=1}^kln(ix+1)\)

根据\(-ln(1-x)=\sum \frac{x^i}{i}\)

原式变为\(\sum_{i=1}^kln(1-(-ix))=\sum_{i=1}^k-\sum_{j=1}^{+∞}\frac{(-ix)^j}{j}\)

即\(-\sum_{i=1}^{+∞}\frac{(-1)^ix^i}{i}\sum_{j=1}^kj^i\)

那么现在需要的就是快速预处理自然数幂和

考虑其\(EGF\)

即\(\sum_{i=1}^{k}\sum_{j=1}^{+∞}\frac{(ix)^j}{j!}=\sum_{i=1}^{k}e^{ix}=\frac{e^{(k+1)x}-e^x}{e^x-1}\)

由于分母常数项为\(0\)，所以上下要同时除以\(x\)

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 N 2000010
#define mo 998244353

vector<int>Is;

int rev[N],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==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(){
	Is.push_back(1);
	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;
	for(int l=1;l<=N-10;l<<=1)G1[l]=mi(3,(mo-1)/(l*2)),G2[l]=mi(G1[l],mo-2);
}

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

struct poly{
	vector<int>val;
	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 modxn(int x){
		if(val.size()<=x)return (poly){val};
		else return (poly){vector<int>(val.begin(),val.begin()+x)};
	}
	poly I(){return (poly){Is};}
	int operator[](int x)const{
		if(x<0||x>=val.size())return 0;
		return val[x];
	}
	void NTT(int x){
		F(i,0,sz()-1)if(rev[i]<i)swap(val[rev[i]],val[i]);
		for(int mid=1;mid<sz();mid<<=1){
			for(int i=0,gn=1;i<sz();i+=(mid<<1),gn=1){
				F(j,0,mid-1){
					int le=val[i+j],ri=1ll*gn*val[i+j+mid]%mo;
					val[i+j]=mod(le+ri);val[i+j+mid]=mod(le-ri+mo);
					gn=1ll*gn*(x==1?G1[mid]:G2[mid])%mo;
				}
			}
		}
		if(x==-1){F(i,0,sz()-1)val[i]=1ll*val[i]*inv[sz()]%mo;}
	}
	void DFT(){NTT(1);}
	void IDFT(){NTT(-1);}
	friend poly operator*(poly x,poly y){
		if(x.sz()<30||y.sz()<30){
			if(x.sz()>y.sz())swap(x,y);
			poly 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);
			}
			return ret;
		}
		int l=1;
		while(l<x.sz()+y.sz()-1)l<<=1;
		x.rsz(l);y.rsz(l);BRT(l);
		x.DFT();y.DFT();
		F(i,0,l-1)x.val[i]=1ll*x[i]*y[i]%mo;
		x.IDFT();
		return x;
	}
	friend poly operator+(poly x,poly y){
		poly 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 operator-(poly x,poly y){
		poly 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 &operator*=(poly x){return (*this)=(*this)*x;}
	poly &operator+=(poly x){return (*this)=(*this)+x;}
	poly &operator-=(poly x){return (*this)=(*this)-x;}
	poly inver(int Len){
		poly f,g;
		f.clear();g.clear();
		g.val.push_back(mi(val[0],mo-2));
		for(int i=2;i<Len*2;i<<=1){
			f.rsz(i<<1);
			g.rsz(i<<1);
			BRT(i<<1);
			F(j,0,i-1)f.val[j]=(j<val.size()?val[j]:0);
			f.DFT();g.DFT();
			F(j,0,(i<<1)-1)g.val[j]=1ll*g[j]*mod(mo+2-1ll*f[j]*g[j]%mo)%mo;
			g.IDFT();
			g.rsz(i);
		}
		return g.modxn(Len);
	}
	poly Ln(int Len){
		poly f,g;
		f.clear();g.clear();
		f.rsz(Len);
		F(i,0,Len-2)f.val[i]=1ll*(i+1>=sz()?0:val[i+1])*(i+1)%mo;
		g=inver(Len);
		g*=f;
		g.modxn(Len);
		Fd(i,Len-1,1)g.val[i]=1ll*g[i-1]*inv[i]%mo;
		g.val[0]=0;
		return g.modxn(Len);
	}
	poly Exp(int Len){
		poly f;
		f.clear();
		f.val.push_back(1);
		for(int i=2;i<Len*2;i<<=1)f=(f*(I()-f.Ln(i)+modxn(i))).modxn(i);
		return f.modxn(Len);
	}
};

int Pk[N],Msum[N];

int main(){
	init();
	int k,m;
	scanf("%d%d",&k,&m);
	poly f,g;
	f.clear();g.clear();
	Pk[0]=1;
	F(i,1,m+1)Pk[i]=1ll*Pk[i-1]*(k+1)%mo;
	F(i,0,m+1)f.val.push_back(1ll*(Pk[i]-1)*ifac[i]%mo),g.val.push_back(ifac[i]);
	F(i,0,m)f.val[i]=f.val[i+1],g.val[i]=g.val[i+1];
	f.val.pop_back();g.val.pop_back();
	g=g.inver(m+1);
	f=(f*g).modxn(m+1);
	F(i,1,m)Msum[i]=1ll*f[i]*fac[i]%mo;
	f.clear();
	f.val.push_back(0);
	F(i,1,m)f.val.push_back(1ll*Msum[i]*((i&1)?1:mo-1)%mo*inv[i]%mo);
	f=f.Exp(m+1);
	F(i,1,m)printf("%d\n",1ll*f[i]*fac[i]%mo);
	return 0;
}

标签：return,加强版,val,int,mo,LuoguP5850,1ll,calc,sum
来源： https://www.cnblogs.com/AmanoKumiko/p/15772412.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

【LuoguP5850】calc加强版