[LOCAL] gcd 求和

题目

题目大意：

对于给定的 (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)) 回答单个询问。

小结：

1. 差分是处理函数（尤其是前缀求和、后缀求和）的有效方式。
2. 思维不应局限，如果反演行不通那么久有必要试试其它方法，例如差分。

代码

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