[LOCAL] gcd 求和

题目

题目大意：

对于给定的 \(n,m\) ，求：

\[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j) \]

数据范围：

\(\text{task_id}\)	\(n,m\le\)	\(T\le\)	特殊性质
\(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\not=m\) ：

  此时可以直接反演，得到结果为：

  \[\sum_{T=1}^{\min\{n,m\}}\lfloor\frac n T\rfloor\lfloor\frac m T\rfloor\varphi(T) \]

• \(n=m\) ：

  此时如果套用反演的结果我们可以愉快地得到：

  \[\sum_{T=1}^n\lfloor\frac n T\rfloor^2\varphi(T) \]

  ......当然是没有办法做的。

  注意到 \(n=m\) 其实是很强的条件，它直接给我们减少了一个变量，因此可以设 \(f(n)=\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)\) 。

  既然直接反演无解，我们不妨用点技巧：对 \(f\) 进行差分。

  \[\begin{aligned} \Delta f(n) &=f(n+1)-f(n)\\ &=2\sum_{i=1}^{n+1}\gcd(i,n+1)-(n+1)\\ &=2\sum_{d|n+1}d\sum_{i=1}^{\frac{n+1}{d}}[(i,d)=1]-(n+1)\\ &=2\sum_{d|n+1}d\varphi(\frac{n+1}{d})-(n+1) \end{aligned} \]

  注意到 \(id*\varphi\) 是积性函数，因此我们可以线筛筛出 \(g(n)=\sum_{d|n}d\varphi(\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( '\n' );
		}
	}
}

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( '\n' );
		}
	}
}

int main()
{
	int task_id;
	read( task_id );
	if( task_id == 6 || task_id == 7 ) Special :: Solve();
	else Normal :: Solve();
	return 0;
}