题目
题目大意：
对于给定的 (n,m) ，求：
[sum_{i=1}^nsum_{j=1}^mgcd(i,j)
]
数据范围：
| ( ext{task_id}) | (n,mle) | (Tle) | 特殊性质 |
|---|---|---|---|
| (1) | (10) | (10^3) | 无 |
| (2) | (10^3) | (10) | 无 |
| (3) | (10^6) | (10^4) | 无 |
| (4) | (10^6) | (10) | (n=m) |
| (5) | (10^6) | (10^4) | (n=m) |
| (6) | (10^6) | (10^5) | (n=m) |
| (7) | (10^7) | (10^6) | (n=m) |
| (8) | (10^6) | (10) | 无 |
| (9) | (10^6) | (10^4) | 无 |
分析
注意数据范围：当 (n=m) 时，数据范围较大，因此我们需要分类计算。
-
(n ot=m) ：
此时可以直接反演，得到结果为：
[sum_{T=1}^{min{n,m}}lfloorfrac n T floorlfloorfrac m T floorvarphi(T) ] -
(n=m) ：
此时如果套用反演的结果我们可以愉快地得到：
[sum_{T=1}^nlfloorfrac n T floor^2varphi(T) ]......当然是没有办法做的。
注意到 (n=m) 其实是很强的条件，它直接给我们减少了一个变量，因此可以设 (f(n)=sum_{i=1}^nsum_{j=1}^ngcd(i,j)) 。
既然直接反演无解，我们不妨用点技巧：对 (f) 进行差分。
[egin{aligned} Delta f(n) &=f(n+1)-f(n)\ &=2sum_{i=1}^{n+1}gcd(i,n+1)-(n+1)\ &=2sum_{d|n+1}dsum_{i=1}^{frac{n+1}{d}}[(i,d)=1]-(n+1)\ &=2sum_{d|n+1}dvarphi(frac{n+1}{d})-(n+1) end{aligned} ]注意到 (id*varphi) 是积性函数，因此我们可以线筛筛出 (g(n)=sum_{d|n}dvarphi(frac n d)) 。因而我们可以 (O(n)) 求出 (f) ，最后 (O(1)) 回答单个询问。
小结：
- 差分是处理函数（尤其是前缀求和、后缀求和）的有效方式。
- 思维不应局限，如果反演行不通那么久有必要试试其它方法，例如差分。
代码
#include <cstdio>
#define rep( i, a, b ) for( int (i) = (a) ; (i) <= (b) ; (i) ++ )
#define per( i, a, b ) for( int (i) = (a) ; (i) >= (b) ; (i) -- )
typedef long long LL;
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;
}
namespace Normal
{
LL val[MAXN];
int pure[MAXN];
int prime[MAXN], pn;
bool isPrime[MAXN];
int N, M;
void EulerSieve( const int siz )
{
val[1] = 1;
for( int i = 2 ; i <= siz ; i ++ )
{
if( ! isPrime[i] )
val[prime[++ pn] = i] = i - 1, pure[i] = 1;
for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= siz ; j ++ )
{
isPrime[i * prime[j]] = true;
if( ! ( i % prime[j] ) )
{
if( ( pure[i * prime[j]] = pure[i] ) == 1 )
val[i * prime[j]] = i * prime[j] - i;
else
val[i * prime[j]] = val[pure[i]] * val[i * prime[j] / pure[i]];
break;
}
val[i * prime[j]] = val[pure[i * prime[j]] = i] * val[prime[j]];
}
}
for( int i = 1 ; i <= siz ; i ++ ) val[i] += val[i - 1];
}
void Solve()
{
int T; EulerSieve( 1e6 );
for( read( T ) ; T -- ; )
{
LL ans = 0;
read( N ), read( M );
for( int l = 1, r ; l <= N && l <= M ; l = r + 1 )
{
r = MIN( N / ( N / l ), M / ( M / l ) );
ans += 1ll * ( N / l ) * ( M / l ) * ( val[r] - val[l - 1] );
}
write( ans ), putchar( '
' );
}
}
}
namespace Special
{
LL val[MAXN], f[MAXN]; int pure[MAXN];
int prime[MAXN], pn;
bool isPrime[MAXN];
int N, M;
void EulerSieve( const int siz )
{
val[1] = 1;
for( int i = 2 ; i <= siz ; i ++ )
{
if( ! isPrime[i] ) val[prime[++ pn] = i] = 2 * i - 1, pure[i] = 1;
for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= siz ; j ++ )
{
isPrime[i * prime[j]] = true;
if( ! ( i % prime[j] ) )
{
if( ( pure[i * prime[j]] = pure[i] ) == 1 )
val[i * prime[j]] = val[i] * prime[j] + 1ll * ( prime[j] - 1 ) * i;
else val[i * prime[j]] = val[pure[i]] * val[i * prime[j] / pure[i]];
break;
}
val[i * prime[j]] = val[pure[i * prime[j]] = i] * val[prime[j]];
}
}
for( int i = 1 ; i <= siz ; i ++ ) f[i] = f[i - 1] + 2 * val[i] - i;
}
void Solve()
{
int T; EulerSieve( 1e7 );
for( read( T ) ; T -- ; )
{
read( N ), read( M );
write( f[N] ), putchar( '
' );
}
}
}
int main()
{
int task_id;
read( task_id );
if( task_id == 6 || task_id == 7 ) Special :: Solve();
else Normal :: Solve();
return 0;
}