题目
分析
简单版
不难发现 (f(n)=mu^2(n)) 。
（这里定义 (f^k(n)=(f(n))^k) ）
然后开始娴熟地推式子：
[egin{aligned}
&sum_{i=1}^nsum_{j=1}^n mu^2(gcd(i,j))gcd(i,j)(i+j)^k\
= &sum_{d=1}^n mu^2(d)dsum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac n d
floor} (id+jd)^k[gcd(i,j)=1]\
= &sum_{d=1}^n mu^2(d)d^{k+1}sum_{t=1}^{lfloorfrac n d
floor}mu(t)t^ksum_{i=1}^{lfloorfrac{n}{dt}
floor}sum_{j=1}^{lfloorfrac{n}{dt}
floor} (i+j)^k\
end{aligned}
]
为了方便简洁，我们定义 (S(n)=sum_{i=1}^nsum_{j=1}^n(i+j)^k) ，然后就有：
[egin{aligned}
&sum_{d=1}^n mu^2(d)d^{k+1}sum_{t=1}^{lfloorfrac n d
floor}mu(t)t^ksum_{i=1}^{lfloorfrac{n}{dt}
floor}sum_{j=1}^{lfloorfrac{n}{dt}
floor} (i+j)^k\
= &sum_{d=1}^n mu^2(d)d^{k+1}sum_{t=1}^{lfloorfrac n d
floor}mu(t)t^kS(lfloorfrac n {dt}
floor)\
end{aligned}
]
现在，只要我们可以快速地求出 (S(n)) ，我们就可以两次整除分块来解决这个问题。
首先不难发现 (S(n)) 有一个相似的定义：
[S(n)=sum_{i=1}^{2n} i^k imes min{i-1,2n-i+1}
]
注意到 (S(n)) 取值只有 (O(sqrt{n})) 种，因此我们可以对于它的每一种取值 (m) ，用 (O(m)) 的时间算出 (S(n)) 。于是我们就可以用 (o(nln n)) 的时间求出所有需要的 (S) 。
然后我们就可以愉快地计算了。
加强版
为了更方便地计算，我们需要继续推式子。下面令 (T=dt) 。
[egin{aligned}
&sum_{d=1}^n mu^2(d)d^{k+1}sum_{t=1}^{lfloorfrac n d
floor}mu(t)S(lfloorfrac n {dt}
floor)\
= &sum_{T=1}^n S(lfloorfrac n T
floor) sum_{d|T} mu^2(d)d^{k+1}mu(frac T d)left(frac{T}{d}
ight)^k\
= &sum_{T=1}^n S(lfloorfrac n T
floor) T^k sum_{d|T} mu^2(d)mu(frac T d)d
end{aligned}
]
式子已经化到底了。现在，如果我们可以快速地求出来 (f(n)=sum_{d|n} mu^2(d)mu(frac n d)d) 和 (S(n)) ，我们就可以运用整除分块，用 (O(sqrt n)) 的时间回答一次查询。
首先考虑 (S) 。不难发现 (S) 实际上是一个三角形系数乘上自然数幂的和（类似于前三角加后三角）。
简单推一推就可以发现，如果令 (g(n)=sum_{i=1}^n (n-i+1) imes i^{k+1}) ，那么 (S(n)=g(2n)-2g(n))
现在我们只需要考虑求 (f(n)) 了。注意到 (f(n)) 实际上是 (mu^2(n)n) 和 (mu(n)) 的狄利克雷卷积，所以它是积性函数。
积性函数，我们就可以使用欧拉筛来求值。考虑新加入一个质因子 (p) ，我们只需要考虑 (f(p^k)) ：
- (f(p^0)=1)
- (f(p^1)=mu^2(1)mu(p)+pmu^2(p)mu(1)=p-1)
- (f(p^2)=mu^2(1)mu(p^2)+pmu^2(p)mu(p)+p^2mu^2(p^2)mu(1)=p)
- (f(p^k)(k>2)=0) ，因为无论怎么分， (mu^2(d)) 和 (mu(frac n d)) 中都会有一个的幂次大于 2 。
于是我们就可以在欧拉筛中，加上判断幂次的过程，并求得 (f) 的值。
注意，本题有那么一点点地卡空间。
一些有价值的点：
- 求出 (f(n)) 的方法比较通用。 (n) 比较小的时候，我们通常可以使用欧拉筛来计算各处的值。（欧拉筛可以被分为增加幂次和增加新质因子两个部分，并且进行分类讨论）
代码
简单版：
#include <cstdio>
typedef long long LL;
const int mod = 998244353;
const int MAXN = 1e7 + 5;
template<typename _T>
void read( _T &x )
{
x = 0;char s = getchar();int f = 1;
while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}
while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}
x *= f;
}
template<typename _T>
void write( _T x )
{
if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }
if( 9 < x ){ write( x / 10 ); }
putchar( x % 10 + '0' );
}
template<typename _T>
_T MIN( const _T a, const _T b )
{
return a < b ? a : b;
}
int su1[MAXN], su2[MAXN];
int pw[MAXN], mu[MAXN];
int prime[MAXN], pn;
bool isPrime[MAXN];
int S[MAXN];
int N; LL K;
int Sub( int x, int v ) { return x < v ? x + mod - v : x - v; }
int Mul( LL x, int v ) { x *= v; if( x >= mod ) x %= mod; return x; }
int Add( int x, int v ) { return x + v >= mod ? x + v - mod : x + v; }
int Qkpow( int base, LL indx )
{
int ret = 1;
while( indx )
{
if( indx & 1 ) ret = Mul( ret, base );
base = Mul( base, base ), indx >>= 1;
}
return ret;
}
void EulerSieve( const int siz )
{
su1[1] = su2[1] = pw[1] = 1;
for( int i = 2 ; i <= siz ; i ++ )
{
if( ! isPrime[i] ) mu[prime[++ pn] = i] = mod - 1, pw[i] = Qkpow( i, K );
for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= siz ; j ++ )
{
isPrime[i * prime[j]] = true, pw[i * prime[j]] = Mul( pw[i], pw[prime[j]] );
if( ! ( i % prime[j] ) ) break; mu[i * prime[j]] = mod - mu[i];
}
su1[i] = Add( su1[i - 1], Mul( Mul( mu[i], mu[i] ), Mul( pw[i], i ) ) );
su2[i] = Add( su2[i - 1], Mul( mu[i], pw[i] ) );
}
}
int GetSum( const int lim )
{
int ret = 0;
for( int k = 1 ; k <= 2 * lim ; k ++ )
ret = Add( ret, Mul( pw[k], MIN( k - 1, 2 * lim - k + 1 ) ) );
return ret;
}
void Init()
{
EulerSieve( N << 1 );
for( int l = 1, r ; l <= N ; l = r + 1 )
{
r = N / ( N / l );
S[N / l] = GetSum( N / l );
}
}
int f( const int n )
{
int ret = 0;
for( int l = 1, r ; l <= n ; l = r + 1 )
{
r = n / ( n / l );
ret = Add( ret, Mul( Sub( su2[r], su2[l - 1] ), S[n / l] ) );
}
return ret;
}
int main()
{
read( N ), read( K );
Init(); int ans = 0;
for( int l = 1, r ; l <= N ; l = r + 1 )
{
r = N / ( N / l );
ans = Add( ans, Mul( Sub( su1[r], su1[l - 1] ), f( N / l ) ) );
}
write( ans ), putchar( '
' );
return 0;
}
加强版：
#include <cstdio>
#include <bitset>
using namespace std;
typedef long long LL;
typedef unsigned int ui;
const int MAXN = 2e7 + 5;
template<typename _T>
void read( _T &x )
{
x = 0;char s = getchar();int f = 1;
while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}
while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}
x *= f;
}
template<typename _T>
void write( _T x )
{
if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }
if( 9 < x ){ write( x / 10 ); }
putchar( x % 10 + '0' );
}
template<typename _T>
_T MIN( const _T a, const _T b )
{
return a < b ? a : b;
}
ui pw[MAXN], f[MAXN];
int prime[MAXN / 10], pn;
bitset<MAXN> isPrime;
int N; LL K;
ui Qkpow( ui base, int indx )
{
ui ret = 1;
while( indx )
{
if( indx & 1 ) ret *= base;
base *= base, indx >>= 1;
}
return ret;
}
void EulerSieve( const int siz )
{
f[1] = pw[1] = 1;
for( int i = 2 ; i <= siz ; i ++ )
{
if( ! isPrime[i] ) pw[prime[++ pn] = i] = Qkpow( i, K ), f[i] = i - 1;
for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= siz ; j ++ )
{
isPrime[i * prime[j]] = true, pw[i * prime[j]] = pw[i] * pw[prime[j]];
if( ! ( i % prime[j] ) )
{
int q = i / prime[j];
if( q % prime[j] ) f[i * prime[j]] = - prime[j] * f[q];
break;
}
f[i * prime[j]] = f[i] * f[prime[j]];
}
}
}
void Init()
{
EulerSieve( N << 1 );
for( int i = 1 ; i <= N << 1 ; i ++ )
f[i] = f[i - 1] + pw[i] * f[i], pw[i] = pw[i - 1] + pw[i];
for( int i = 1 ; i <= N << 1 ; i ++ ) pw[i] = pw[i - 1] + pw[i];
}
int S( const int n ) { return pw[n << 1] - ( pw[n] << 1 ); }
int main()
{
int T, n;
read( T ), read( N ), read( K );
Init();
while( T -- )
{
read( n ); ui ans = 0;
for( int l = 1, r ; l <= n ; l = r + 1 )
{
r = n / ( n / l );
ans += ( f[r] - f[l - 1] ) * S( n / l );
}
write( ans ), putchar( '
' );
}
return 0;
}