Solution -「洛谷 P6156」简单题

Description

Link.

求 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^n(i+j)^kf(\gcd(i,j))\gcd(i,j)\)。

Solution

\[\begin{aligned} \textbf{ANS}&=\sum_{i=1}^{n}\sum_{j=1}^{n}(i+j)^{k}\mu^{2}(\gcd(i,j))\gcd(i,j) \\ &=\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{n}(i+j)^{k}\mu^{2}(d)d[\gcd(i,j)=d] \\ &=\sum_{d=1}^{n}d^{k+1}\times\mu^{2}(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}(i+j)^{k}\sum_{h|i,h|j}\mu(h) \\ &=\sum_{d=1}^{n}d^{k+1}\times\mu^{2}(d)\sum_{h=1}^{\lfloor\frac{n}{d}\rfloor}\mu(h)\times h^{k}\times\sum_{i=1}^{\lfloor\frac{n}{dh}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dh}\rfloor}(i+j)^{k} \\ &=\sum_{d=1}^{n}d^{k+1}\times\mu^{2}(d)\sum_{h=1}^{\lfloor\frac{n}{d}\rfloor}\mu(h)\times h^{k}\times\sum_{i=1}^{\lfloor\frac{n}{dh}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dh}\rfloor}(i+j)^{k} \\ \end{aligned} \\ \]

前面两个和式里面显然能算，考虑怎么对于 \(x\) 算 \(\sum_{i=1}^{x}\sum_{j=1}^{x}(i+j)^{k}\)。考虑对其差分：

\[\begin{aligned} \left(\sum_{i=1}^{x+1}\sum_{j=1}^{x+1}(i+j)^{k}\right)-\left(\sum_{i=1}^{x}\sum_{j=1}^{x}(i+j)^{k}\right)&=\sum_{i=1}^{x}\sum_{j=1}^{x+1}(i+j)^{k}+\sum_{i=1}^{x+1}(x+1+i)^{k}-\sum_{i=1}^{x}\sum_{j=1}^{x}(i+j)^{k} \\ &=\sum_{i=1}^{x}\left(\sum_{j=1}^{x+1}(i+j)^{k}-\sum_{j=1}^{x}(i+j)^{k}\right)+\sum_{i=1}^{x+1}(x+1+i)^{k} \\ &=\sum_{i=1}^{x}(x+1+i)^{k}+\sum_{i=1}^{x+1}(x+1+i)^{k} \\ \end{aligned} \]

然后滚个前缀和就可以算了。

#include<bits/stdc++.h>
typedef long long LL;
const int MOD=998244353;
int norm( LL x ) {
	if( x<0 ) {
		x+=MOD;
	}
	if( x>=MOD ) {
		x%=MOD;
	}
	return x;
}
int n,k,ans;
int qpow( int bas,int fur ) {
	int res=1;
	while( fur ) {
		if( fur&1 ) {
			res=norm( LL( res )*bas );
		}
		bas=norm( LL( bas )*bas );
		fur>>=1;
	}
	return norm( res+MOD );
}
std::tuple<std::vector<int>,std::vector<int>> makePrime( int n ) {
	std::vector<int> prime,tag( n+1 ),mu( n+1 ),pw( n+1 );
	pw[0]=1;
	mu[1]=pw[1]=1;
	for( int i=2;i<=n;++i ) {
		if( !tag[i] ) {
			mu[i]=norm( -1 );
			prime.emplace_back( i );
			pw[i]=qpow( i,k );
		}
		for( int j=0;j<int( prime.size() ) && i*prime[j]<=n;++j ) {
			tag[i*prime[j]]=1;
			pw[i*prime[j]]=norm( LL( pw[i] )*pw[prime[j]] );
			if( i%prime[j]==0 ) {
				mu[i*prime[j]]=0;
				break;
			} else {
				mu[i*prime[j]]=norm( -mu[i] );
			}
		}
	}
	return std::tie( mu,pw );
}
int main() {
	LL tmp;
	scanf( "%d %lld",&n,&tmp );
	k=tmp%( MOD-1 );
	std::vector<int> mu,pw,prt( n+1 ),exprt( n+1 ),preSum( n+1 );
	// prt: i^(k+1)*mu^2(i)
	// exprt: mu(i)*i^k
	// preSum sum sum (i+j)^k
	std::tie( mu,pw )=makePrime( n<<1|1 );
	for( int i=1;i<=n;++i ) {
		prt[i]=norm( prt[i-1]+norm( LL( norm( LL( norm( LL( mu[i] )*mu[i] ) )*pw[i] ) )*i ) );
		exprt[i]=norm( exprt[i-1]+norm( LL( mu[i] )*pw[i] ) );
	}
	for( int i=1;i<=( n<<1 );++i ) {
		pw[i]=norm( pw[i]+pw[i-1] );
	}
	for( int i=1;i<=n;++i ) {
		preSum[i]=norm( norm( preSum[i-1]+norm( pw[i<<1]-pw[i] ) )+norm( pw[(i<<1)-1]-pw[i] ) );
	}
	for( int l=1,r;l<=n;l=r+1 ) {
		r=n/( n/l );
		int tmp=0;
		for( int exl=1,exr,m=n/l;exl<=m;exl=exr+1 ) {
			exr=m/( m/exl );
			tmp=norm( tmp+norm( LL( norm( exprt[exr]-exprt[exl-1] ) )*preSum[m/exl] ) );
		}
		ans=norm( ans+LL( norm( prt[r]-prt[l-1] ) )*tmp );
	}
	printf("%d\n",ans);
	return 0;
}