题解 P5850 【calc加强版】

题意

一个序列\(a_1,a_2,\ldots,a_n\)合法，当且仅当：

• \(\forall i\in[1,n],a_i\in[1,k]\)
• \(\forall i\not=j,a_i\not=a_j\)

一个序列的值为它的乘积,即\(\prod_{i=1}^na_i\)

给定\(m,k\),对于\(n\in[1,m]\)，求所有长度为\(n\)的不同的序列的值的和。

题解

不难发现答案可以写成生成函数的形式：

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

但此时是无序的，因为\(a_1,a_2,\ldots,a_n\)互不相等，乘\(n!\)即可。最后答案为\(n![x^n]F(x)\)

可以写出\(\ln F(x)\)，证明见上一篇题解\(\color{skyblue}{\text{P4463 【[集训队互测2012] calc】}}\)

\[\ln F(x)=\sum_{j=1}^\infty\frac{(-1)^{j-1}}{j}(\sum_{i=1}^ki^j)x^j \]

所以我们现在希望能对于\(j\in[1,n]\)求出\(\sum_{i=1}^ki^j\)的值。

考虑\(<\sum_{i=1}^ki^0,\sum_{i=1}^ki^1,\sum_{i=1}^ki^2,\ldots>\)的\(\mathbf{EFG}\)

\[G(x)=\sum_{j=0}^\infty\frac{\sum_{i=1}^ki^j}{j!}x^j \]

我们交换一下求和顺序

\[=\sum_{i=1}^k\sum_{j=0}^\infty \frac{(ix)^j}{j!} \]

右边那串实际上就是\(e^{ix}\)的泰勒展开。

\[=\sum_{i=1}^ke^{ix} \]

发现就是一个长度为\(k\)，首项为\(e^x\)公比为\(e^x\)的等差数列，套用公式得到：

\[=\frac{e^{(k+1)x}-e^x}{e^x-1} \]

但此时分子的常数项项为\(0\)不存在乘法逆。\(\rm However,\)分母的常数项也为\(0\)因此上下同除\(x\)就可以放心大胆地求乘法逆了。

代码

#include<bits/stdc++.h>
namespace in{
	char buf[1<<21],*p1=buf,*p2=buf;
	inline int getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
	template <typename T>inline void read(T& t){
		t=0;int f=0;char ch=getc();while (!isdigit(ch)){if(ch=='-')f = 1;ch=getc();}
	    while(isdigit(ch)){t=t*10+ch-48;ch = getc();}if(f)t=-t;
	}
	template <typename T,typename... Args> inline void read(T& t, Args&... args){read(t);read(args...);}
}
namespace out{
	char buffer[1<<21];int p1=-1;const int p2 = (1<<21)-1;
	inline void flush(){fwrite(buffer,1,p1+1,stdout),p1=-1;}
	inline void putc(const char &x) {if(p1==p2)flush();buffer[++p1]=x;}
	template <typename T>void write(T x) {
		static char buf[15];static int len=-1;if(x>=0){do{buf[++len]=x%10+48,x/=10;}while (x);}else{putc('-');do {buf[++len]=-(x%10)+48,x/=10;}while(x);}
		while (len>=0)putc(buf[len]),--len;
	}
}
using namespace std;
template<const int mod>
struct modint{
    int x;
    modint<mod>(int o=0){x=o;}
    modint<mod> &operator = (int o){return x=o,*this;}
    modint<mod> &operator +=(modint<mod> o){return x=x+o.x>=mod?x+o.x-mod:x+o.x,*this;}
    modint<mod> &operator -=(modint<mod> o){return x=x-o.x<0?x-o.x+mod:x-o.x,*this;}
    modint<mod> &operator *=(modint<mod> o){return x=1ll*x*o.x%mod,*this;}
    modint<mod> &operator ^=(int b){
        modint<mod> a=*this,c=1;
        for(;b;b>>=1,a*=a)if(b&1)c*=a;
        return x=c.x,*this;
    }
    modint<mod> &operator /=(modint<mod> o){return *this *=o^=mod-2;}
    modint<mod> &operator +=(int o){return x=x+o>=mod?x+o-mod:x+o,*this;}
    modint<mod> &operator -=(int o){return x=x-o<0?x-o+mod:x-o,*this;}
    modint<mod> &operator *=(int o){return x=1ll*x*o%mod,*this;}
    modint<mod> &operator /=(int o){return *this *= ((modint<mod>(o))^=mod-2);}
	template<class I>friend modint<mod> operator +(modint<mod> a,I b){return a+=b;}
    template<class I>friend modint<mod> operator -(modint<mod> a,I b){return a-=b;}
    template<class I>friend modint<mod> operator *(modint<mod> a,I b){return a*=b;}
    template<class I>friend modint<mod> operator /(modint<mod> a,I b){return a/=b;}
    friend modint<mod> operator ^(modint<mod> a,int b){return a^=b;}
    friend bool operator ==(modint<mod> a,int b){return a.x==b;}
    friend bool operator !=(modint<mod> a,int b){return a.x!=b;}
    bool operator ! () {return !x;}
    modint<mod> operator - () {return x?mod-x:0;}
	modint<mod> &operator++(int){return *this+=1;}
};
const int N=4e6+5;

const int mod=998244353;
namespace residue{
	//求二次剩余
	int I2;
	struct complex{
		long long real,imag;
		complex&operator=(long long x){real=x,imag=0;}
		complex(long long real=0,long long imag=0):real(real),imag(imag){}
		bool operator==(const complex&b)const{return real==b.real&&imag==b.imag;}
		complex operator*(const complex &b)const{return complex((real*b.real+I2*imag%mod*b.imag)%mod,(real*b.imag+imag*b.real)%mod);}
	};
	complex ksm(complex a,int b){
		complex res=1;
		while(b){
			if(b&1)res=res*a;
			a=a*a;b>>=1;
		}return res;
	}
	bool check(int x){return ksm(complex(x,0),(mod-1)/2).real==1;}
	pair<int,int>solve(int n){
	    if(n==0)return {0,0};
		long long a=rand()%mod;
		while(!a||check((a*a+mod-n)%mod))a=rand()%mod;
		I2=(a*a+mod-n)%mod;
		int x0=ksm(complex(a,1),(mod+1)/2).real;
		return {min(x0,mod-x0),max(x0,mod-x0)};
	}
}
const modint<mod> GG=3,Ginv=modint<mod>(1)/3,I=86583718;
struct poly{
	vector<modint<mod>>a;
	modint<mod>&operator[](int i){return a[i];}
	int size(){return a.size();}
	void resize(int n){a.resize(n);}
	void reverse(){std::reverse(a.begin(),a.end());}
};
int rev[N];
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){int n=1<<k;for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));}
inline void ntt(poly&a,int k,int typ){
	int n=1<<k;
	for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int mid=1;mid<n;mid<<=1){
		modint<mod> wn=(typ>0?GG:Ginv)^((mod-1)/(mid<<1));
		for(int r=mid<<1,j=0;j<n;j+=r){
			modint<mod> w=1;
			for(int k=0;k<mid;k++,w=w*wn){
				modint<mod> x=a[j+k],y=w*a[j+k+mid];
				a[j+k]=x+y,a[j+k+mid]=x-y;
			}
		}
	}
	if(typ<0){
		modint<mod> inv=modint<mod>(1)/n;
		for(int i=0;i<n-1;i++)a[i]*=inv;
	}
}
inline poly one(){poly a;a.a.push_back(1);return a;}
poly operator +(poly a,poly b){
	int n=max(a.size(),b.size());a.resize(n),b.resize(n);
	for(int i=0;i<n;i++)a[i]+=b[i];return a;
}
poly operator -(poly a,poly b){
	int n=max(a.size(),b.size());a.resize(n),b.resize(n);
	for(int i=0;i<n;i++)a[i]-=b[i];return a;
}
inline poly operator*(poly a,poly b){
	int n=a.size()+b.size()-1,k=ext(n);
	a.resize(1<<k),b.resize(1<<k),init(k);
	ntt(a,k,1);ntt(b,k,1);for(int i=0;i<(1<<k);i++)a[i]*=b[i];
	ntt(a,k,-1),a.resize(n);return a;
}
inline poly operator*(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]*=b;return a; }
inline poly operator/(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]/=b;return a; }
inline poly operator-(poly a){for(int i=0;i<a.size();i++)a[i]=-a[i];return a; }
poly inv(poly F,int k){
	int n=1<<k;F.resize(n);
	if(n==1){F[0]=modint<mod>(1)/F[0];return F;}
	poly G,H=inv(F,k-1);
	G.resize(n),H.resize(n<<1),F.resize(n<<1);
	for(int i=0;i<n/2;i++)G[i]=H[i]*2;
	init(k+1),ntt(H,k+1,1),ntt(F,k+1,1);
	for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i]*F[i];
	ntt(H,k+1,-1),H.resize(n);
	for(int i=0;i<n;i++)G[i]-=H[i];return G;
}
inline poly inv(poly a){
	int n=a.size();
	a=inv(a,ext(n)),a.resize(n);return a;;
}
inline poly deriv(poly a){//求导 
	int n=a.size()-1;
	for(int i=0;i<n;i++)a[i]=a[i+1]*(i+1);
	a.resize(n);return a;
}
inline poly inter(poly a){//求原 
	int n=a.size()+1;a.resize(n);
	for(int i=n;i>=1;i--)a[i]=a[i-1]/i;
	a[0]=0;return a;
}
inline poly ln(poly a){
	int n=a.size();
	a=inter(deriv(a)*inv(a));
	a.resize(n);return a;
}
poly exp(poly a,int k){
	int n=1<<k;a.resize(n);
	if(n==1)return one();
	poly f0=exp(a,k-1);f0.resize(n);
	return f0*(one()+a-ln(f0)); 
}
poly exp(poly a){
	int n=a.size();
	a=exp(a,ext(n));a.resize(n);return a;
}
typedef modint<mod>mint;
mint fac[N];
void add(poly &F,mint a,bool flag){
	//flag?F+=e^ax:F-=e^ax
	for(int i=0;i<F.size();i++)
		flag?(F[i]+=(a^i)/fac[i]):(F[i]-=(a^i)/fac[i]);
}
int n,k;poly F,G;
signed main(){
	fac[0]=1;for(int i=1;i<N;i++)fac[i]=fac[i-1]*i;
	in::read(k,n);
	F.resize(n+2);G.resize(n+2); 
	add(F,k+1,1);add(F,1,0);add(G,1,1);
	for(int i=0;i<=n;i++)F[i]=F[i+1],G[i]=G[i+1];
	F.resize(n+1);G.resize(n+1); 
	F=F*inv(G);G[0]=0;
	for(int i=1;i<=n;i++)
		if((i+1)&1)G[i]=-F[i]*fac[i]/i;
		else G[i]=F[i]*fac[i]/i;
	G=exp(G);for(int i=1;i<=n;i++)out::write((G[i]*fac[i]).x),out::putc('\n');
	out::flush();
	return 0;
}