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;
}