Description
Link.
求 (sumlimits_{i=1}^nsumlimits_{j=1}^n(i+j)^kf(gcd(i,j))gcd(i,j))。
Solution
[egin{aligned}
extbf{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} imesmu^{2}(d)sum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{n}{d}
floor}(i+j)^{k}sum_{h|i,h|j}mu(h) \
&=sum_{d=1}^{n}d^{k+1} imesmu^{2}(d)sum_{h=1}^{lfloorfrac{n}{d}
floor}mu(h) imes h^{k} imessum_{i=1}^{lfloorfrac{n}{dh}
floor}sum_{j=1}^{lfloorfrac{n}{dh}
floor}(i+j)^{k} \
&=sum_{d=1}^{n}d^{k+1} imesmu^{2}(d)sum_{h=1}^{lfloorfrac{n}{d}
floor}mu(h) imes h^{k} imessum_{i=1}^{lfloorfrac{n}{dh}
floor}sum_{j=1}^{lfloorfrac{n}{dh}
floor}(i+j)^{k} \
end{aligned} \
]
前面两个和式里面显然能算，考虑怎么对于 (x) 算 (sum_{i=1}^{x}sum_{j=1}^{x}(i+j)^{k})。考虑对其差分：
[egin{aligned}
left(sum_{i=1}^{x+1}sum_{j=1}^{x+1}(i+j)^{k}
ight)-left(sum_{i=1}^{x}sum_{j=1}^{x}(i+j)^{k}
ight)&=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}
ight)+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
",ans);
return 0;
}