[洛谷 P6156] 简单题 及其加强
题目
分析
简单版
不难发现 \(f(n)=\mu^2(n)\) 。
（这里定义 \(f^k(n)=(f(n))^k\) ）
然后开始娴熟地推式子：
\[\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^n \mu^2(\gcd(i,j))\gcd(i,j)(i+j)^k\\
= &\sum_{d=1}^n \mu^2(d)d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor} (id+jd)^k[\gcd(i,j)=1]\\
= &\sum_{d=1}^n \mu^2(d)d^{k+1}\sum_{t=1}^{\lfloor\frac n d\rfloor}\mu(t)t^k\sum_{i=1}^{\lfloor\frac{n}{dt}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dt}\rfloor} (i+j)^k\\
\end{aligned}
\]
为了方便简洁，我们定义 \(S(n)=\sum_{i=1}^n\sum_{j=1}^n(i+j)^k\) ，然后就有：
\[\begin{aligned}
&\sum_{d=1}^n \mu^2(d)d^{k+1}\sum_{t=1}^{\lfloor\frac n d\rfloor}\mu(t)t^k\sum_{i=1}^{\lfloor\frac{n}{dt}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dt}\rfloor} (i+j)^k\\
= &\sum_{d=1}^n \mu^2(d)d^{k+1}\sum_{t=1}^{\lfloor\frac n d\rfloor}\mu(t)t^kS(\lfloor\frac n {dt}\rfloor)\\
\end{aligned}
\]
现在，只要我们可以快速地求出 \(S(n)\) ，我们就可以两次整除分块来解决这个问题。
首先不难发现 \(S(n)\) 有一个相似的定义：
\[S(n)=\sum_{i=1}^{2n} i^k\times \min\{i-1,2n-i+1\}
\]
注意到 \(S(n)\) 取值只有 \(O(\sqrt{n})\) 种，因此我们可以对于它的每一种取值 \(m\) ，用 \(O(m)\) 的时间算出 \(S(n)\) 。于是我们就可以用 \(o(n\ln n)\) 的时间求出所有需要的 \(S\) 。
然后我们就可以愉快地计算了。
加强版
为了更方便地计算，我们需要继续推式子。下面令 \(T=dt\) 。
\[\begin{aligned}
&\sum_{d=1}^n \mu^2(d)d^{k+1}\sum_{t=1}^{\lfloor\frac n d\rfloor}\mu(t)S(\lfloor\frac n {dt}\rfloor)\\
= &\sum_{T=1}^n S(\lfloor\frac n T\rfloor) \sum_{d|T} \mu^2(d)d^{k+1}\mu(\frac T d)\left(\frac{T}{d}\right)^k\\
= &\sum_{T=1}^n S(\lfloor\frac n T\rfloor) 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)\times 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)+p\mu^2(p)\mu(1)=p-1\)
- \(f(p^2)=\mu^2(1)\mu(p^2)+p\mu^2(p)\mu(p)+p^2\mu^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( '\n' );
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( '\n' );
}
return 0;
}
