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^jsum_{i=1}^{frac{n}{d}}i^j ]

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

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

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

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

暴力带进去。

[=sum_{d|n}mu(d)sum_{j=0}^{m-1}a_jd^jsum_{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^jsum_{i=1}^{frac{n}{d}}i^j\ ]

多算 (0^0) 当且仅当 (j=0) ,这时候贡献是 (sum_{d|n}mu(d)a_0=a_0 imes[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^jsum_{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
",B[i]);
	return 0;
}
原文地址:https://www.cnblogs.com/zzctommy/p/14301079.html