P7102 [w3R1] 算

P7102 [w3R1] 算

\[s(n)=\sum_{i=1}^{n}p(i)\sum_{d|\gcd(i,n)}\mu(d)\\ =\sum_{d|n}\mu(d)\sum_{i=1}^{\frac{n}{d}}p(id)\\ =\sum_{d|n}\mu(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=0}^{m-1}a_j(id)^{j}\\ =\sum_{d|n}\mu(d)\sum_{j=0}^{m-1}a_jd^j\sum_{i=1}^{\frac{n}{d}}i^j \]

后面有一个自然数幂和的东西，考虑直接伯努利数求通项，这是一个 \(j+1\) 次多项式。设这个多项式为

\[S_k(n)=\sum\limits_{i=0}^{n}i^k=k!\sum\limits_{i=1}^{k+1}\dfrac{B_i}{(k+1-i)!}n^{k+1-i}+n^k \]

如果你不了解怎么搞这个通项，可以看 多项式笔记（二） 。（这里的 \(B_i\) 是伯努利数 \(\rm EGF\) 的 \([x^i]\) 系数，也就是除以了阶乘的）

不过这里注意一件事情，伯努利数的求和上界是 \(n-1\) ，所以用伯努利数算 \(B\) 的时候，最后要单独加 \(n^k\) 。以及伯努利数默认 \(0^0=1\) ，要减掉 。这个一般弱于前面的那个和式，单独讨论一下就能解决，留到最后再说。

暴力带进去。

\[=\sum_{d|n}\mu(d)\sum_{j=0}^{m-1}a_jd^j\sum_{i=1}^{\frac{n}{t}}i^j\\ =\sum_{d|n}\mu(d)\sum_{j=0}^{m-1}a_jd^jj!\sum_{i=1}^{j+1}\dfrac{B_i}{(j+1-i)!}(\dfrac{n}{d})^{j+1-i}\\ =\sum_{j=0}^{m-1}a_jj!\sum_{i=0}^{j}\dfrac{B_i}{(j+1-i)!}n^{j+1-i}\sum_{d|n}\mu(d)d^{i-1}\\ \]

后面那个是积性函数，看几眼就感觉很可做。

考虑把 \(n\) 质因数分解了，对于每一个 \(p^k\) 算答案，然后乘起来。

观察到 \(\mu(p^k)\) 当且仅当在 \(k=0\) 或 \(k=1\) 时有值。

所以

\[\sum_{d|p^k}\mu(d)d^{i-1}=1-p*p^{i-1}=1-p^i \]

令 \(f(k,n)=\sum_{d|n}\mu(d)d^k\) ，那么显然有 \(f(k,n^c)=f(k,n)\) ，因为一旦存在 \(p\) 这个因子就必然只会产生 \(1-p^k\) 的贡献，与因子数量无关。

Pollard-Rho预处理质因数，可以 \(O(\log n)\) 左右算出来，不是瓶颈就是了。

那么可以直接预处理所有 \(A_i=f(i-1,c)\)

\[=\sum_{j=0}^{m-1}a_jj!\sum_{i=0}^{j}n^{j+1-i}\dfrac{B_{i}}{(j+1-i)!}A_i\\ =\sum_{j=0}^{m-1}a_jj!\sum_{i=0}^{j}n^{i+1}\dfrac{B_{j-i+1}}{(i+1)!}A_{j-i+1}\\ =\sum_{i=0}^{m-1}n^{i+1}\dfrac{1}{(i+1)!}\sum_{j=i}^{m-1}a_{j}j!A_{j-i}B_{j-i} \]

发现这个就是一个与 \(n\) 有关的多项式，系数差卷积一下就可以出来。然后我们需要求出连续的 \(c^i\) 的点值，就是CZT变换板子，不会可以看 多项式笔记（一）。

然后来讨论伯努利数边界的问题：我们会多算所有 \(0^0\) ，要减掉；会漏算所有 \((\dfrac{n}{t})^j\) ，要加上。

拉一个比较原始的式子下来：

\[\sum_{d|n}\mu(d)\sum_{j=0}^{m-1}a_jd^j\sum_{i=1}^{\frac{n}{d}}i^j\\ \]

多算 \(0^0\) 当且仅当 \(j=0\) ，这时候贡献是 \(\sum_{d|n}\mu(d)a_0=a_0\times[n=1]\) ，也就是说只有在 \(n=1\) 的时候会多算 \(a_0\) 。

这个好办，直接判一下 \(c=1\) 的情况，输出 \(t\) 个 \(p(1)\) 即可。

少算 \((\dfrac{n}{d})^j\) 稍微棘手一些

\[\sum_{d|n}\mu(d)\sum_{j=0}^{m-1}a_jd^j(\dfrac{n}{d})^j\\ =\sum_{d|n}\mu(d)\sum_{j=0}^{m-1}a_jn^j\\ =\sum_{j=0}^{m-1}a_jn^j\sum_{d|n}\mu(d)\\ =\sum_{j=0}^{m-1}a_jn^j[n=1]\\ \]

发现还是在 \(n=1\) 的时候会出事，然而已经判掉了所以啥事没有。

好难码啊，边界乱七八糟的/kk。

不过出题人貌似不想给完整代码（详见commd_block洛谷题解评论），所以我也不给完整代码了/cy。但是删去部分都是很简单的部分，只防xxs，并且加了详细注释，如果你有做这题的水平必然看得懂。

LL mul(LL x,LL y,LL mod){
	LL res=x*y-(long long)((long double)x/mod*y+0.5)*mod;
	return res<0?res+mod:res;
}
LL gcd(LL x,LL y){
	while(y){LL t=x%y;x=y,y=t;}
	return x;
}
LL qpow(LL n,LL k,LL mod){
	LL res=1;for(;k;k>>=1,n=mul(n,n,mod))if(k&1)res=mul(res,n,mod);
	return res;
}

namespace MR{
int testp[8]={2,3,5,7,11,19,61,233};
bool mr(LL x,LL p){
	if(x%p==0||qpow(p,x-1,x)!=1)return 0;
	LL k=x-1;
	while(!(k&1)){
		LL t=qpow(p,k>>=1,x);
		if(t!=x-1&&t!=1)return 0;
		if(t==x-1)break;
	}
	return 1;
}
bool test(LL x){
	for(int i=0;i<8;++i){
		if(x==testp[i])return 1;
		if(!mr(x,testp[i]))return 0;
	}
	return 1;
}
}

namespace PR{
int tot;
LL d[150];
mt19937_64 rnd(673);
LL pr(LL x,LL c){
	if(x==4)return 2;
	static LL s[150],v0,v1,g;
	static int len;
	len=0,v0=c,v1=mul(v0,v0,x)+c,g=1;
	while(1){
		s[++len]=v1-v0,g=mul(g,v1-v0,x);
		if(len==127){if(gcd(g,x)>1)break;len=0;}
		v0=mul(v0,v0,x)+c,v1=mul(v1,v1,x)+c,v1=mul(v1,v1,x)+c;
	}
	for(int i=1;i<=len;++i)if((g=gcd(s[i],x))>1)return g;
	return x;
}
void solve(LL x){
	if(x==1)return;
	if(MR::test(x))return d[++tot]=x,void();
	LL y=x;while(y==x)y=pr(x,rnd()%(x-1)+1);
	solve(x/y),solve(y);
}
void work(LL x){tot=0,solve(x),sort(d+1,d+tot+1);}

}

const int N=200005;
const int M=N*3*2;
#define mod 998244353

namespace math{
int inv[N],fac[N],ifc[N];
inline int qpow(int n,int k){int res=1;for(;k;k>>=1,n=1ll*n*n%mod)if(k&1)res=1ll*n*res%mod;return res;}
inline int comb(int n,int m){return n<m?0:1ll*fac[n]*ifc[m]%mod*ifc[n-m]%mod;}
inline void fmod(int&x){x-=mod,x+=x>>31&mod;}
void initmath(const int&n=N-1){
	fac[0]=1;for(int i=1;i<=n;++i)fac[i]=1ll*fac[i-1]*i%mod;
	ifc[n]=qpow(fac[n],mod-2);for(int i=n-1;i>=0;--i)ifc[i]=1ll*ifc[i+1]*(i+1)%mod;
	inv[1]=1;for(int i=2;i<=n;++i)inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
}
using namespace math;

namespace poly{

int rev[M],rt[M],lim;
void initpoly(const int&n){
	static int lg;
	for(lg=0,lim=1;lim<n;++lg,lim<<=1);
	for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
	for(int i=1;i<lim;i<<=1){
		int w=qpow(3,(mod-1)/(i<<1));
		rt[i]=1;
		for(int j=1;j<i;++j)rt[i+j]=1ll*rt[i+j-1]*w%mod;
	}
}
void NTT(int*a,int op){
	if(!op)reverse(a+1,a+lim);
	for(int i=0;i<lim;++i)
		if(i>rev[i])swap(a[i],a[rev[i]]);
	for(int i=1;i<lim;i<<=1){
		for(int j=0;j<lim;j+=i<<1){
			for(int k=0;k<i;++k){
				const int t=1ll*rt[i+k]*a[i+j+k]%mod;
				fmod(a[i+j+k]=a[j+k]+mod-t),fmod(a[j+k]+=t);
			}
		}
	}
	if(!op)for(int i=0,iv=qpow(lim,mod-2);i<lim;++i)a[i]=1ll*a[i]*iv%mod;
}
#define clr(a,n) memset(a,0,sizeof(int)*(n))
#define cpy(a,b,n) memcpy(a,b,sizeof(int)*(n))
void poly_mul(int*f,int*g,int*ans,int n,int m){//这是卷积，f[0,1,...,n-1]*g[0,1,...,m-1]->ans[0,1,...,n+m-2]
}
void poly_inv(int*g,int*f,int n){//这是多项式求逆，f[0,1,...,n-1]的逆元->g[0,1,...,n-1]
}
void bernoulli(int*g,int n){
	static int A[M];
	for(int i=0;i<n;++i)A[i]=ifc[i+1];
	clr(g,n),poly_inv(g,A,n);
}
void CZT(int*g,int*f,int m,int n,int c){//给定n-1次多项式f，求f(c^{0->m-1})
	static int pw[N<<1],ipw[N<<1],ivc,A[M],B[M];
	pw[0]=pw[1]=ipw[0]=ipw[1]=1,ivc=qpow(c,mod-2);
	for(int i=2;i<n+m;++i)pw[i]=1ll*pw[i-1]*c%mod,ipw[i]=1ll*ipw[i-1]*ivc%mod;
	for(int i=2;i<n+m;++i)pw[i]=1ll*pw[i-1]*pw[i]%mod,ipw[i]=1ll*ipw[i-1]*ipw[i]%mod;
	for(int i=0;i<n;++i)A[i]=1ll*f[i]*ipw[i]%mod;
	for(int i=0;i<n+m;++i)B[i]=pw[i];
	reverse(A,A+n),poly_mul(A,B,A,n,n+m);
	for(int i=0;i<m;++i)g[i]=1ll*ipw[i]*A[i+n-1]%mod;
}

}

using PR::d;
using PR::tot;
int m,t,a[N],A[M],B[M],tem[N],num,pri[N];
LL c;
signed main(){
	initmath();
	scanf("%d%lld%d",&m,&c,&t),++t;
	for(int i=0;i<m;++i)a[i]=read();
	if(c==1){
		static LL sum;
		sum=0;for(int i=0;i<m;++i)sum+=a[i];
		int s=sum%mod;
		for(int i=1;i<t;++i)printf("%d ",s);
		return 0;
	}
	PR::work(c);
	for(int l=1,r=0;l<=tot;l=r+1){
		while(r<tot&&d[r+1]==d[l])++r;
		pri[++num]=d[l]%mod;
	}
	A[0]=1;
	for(int i=1;i<=num;++i)A[0]=1ll*A[0]*(mod+1-(tem[i]=qpow(pri[i],mod-2)))%mod;
	for(int i=1;i<m;++i){
		A[i]=1;
		for(int j=1;j<=num;++j)A[i]=1ll*A[i]*(mod+1-(tem[j]=1ll*tem[j]*pri[j]%mod))%mod;
	}
	poly::bernoulli(B,m);
	for(int i=0;i<m;++i)A[i]=1ll*A[i]*B[i]%mod;
	for(int i=0;i<m;++i)B[i]=1ll*a[i]*fac[i]%mod;
	reverse(A,A+m),poly::poly_mul(A,B,B,m,m);
	for(int i=0;i<m;++i)A[i+1]=1ll*ifc[i+1]*B[i+m-1]%mod;
	A[0]=0,clr(B,t),poly::CZT(B,A,t,m+1,c%mod);
	for(int i=1;i<t;++i)printf("%d\n",B[i]);
	return 0;
}