Solution -「洛谷 P4449」于神之怒加强版

\(\mathcal{Description}\)

  Link.

  给定 \(k\) 和 \(T\) 组 \(n,m\)，对于每组，求

\[\sum_{i=1}^n\sum_{j=1}^m\operatorname{gcd}^k(i,j)\bmod(10^9+7) \]

  \(T\le2\times10^3\)，\(n,m,k\le5\times10^6\)。

\(\mathcal{Solution}\)

  几个月没推式子找找手感 qwq。（

  不妨设 \(n\le m\)：

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^m\operatorname{gcd}^k(i,j)&=\sum_{d=1}^nd^k\sum_{i=1}^{\lfloor\frac{n}d \rfloor}\sum_{j=1}^{\lfloor\frac{m}d \rfloor}[i\perp j]\\ &=\sum_{d=1}^nd^k\sum_{d'=1}^{\lfloor\frac{n}d \rfloor}\mu(d')\lfloor\frac{n}{dd'}\rfloor\lfloor\frac{m}{dd'}\rfloor\\ &=\sum_{T=1}^n\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d|T}d^k\mu(\frac{T}d),~~~~\text{let }T=dd'\\ &=\sum_{T=1}^n\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor(\operatorname{id}^k\star\mu)(T) \end{aligned} \]

  \(\operatorname{id}^k\star\mu\) 积性，可以线性筛筛出。此后整除分块处理询问。复杂度 \(\mathcal O(n)-\mathcal O(\sqrt n)\)。

\(\mathcal{Code}\)

/* Clearink */

#include <cstdio>

#define rep( i, l, r ) for ( int i = l, repEnd##i = r; i <= repEnd##i; ++i )
#define per( i, r, l ) for ( int i = r, repEnd##i = l; i >= repEnd##i; --i )

inline int rint () {
	int x = 0, f = 1; char s = getchar ();
	for ( ; s < '0' || '9' < s; s = getchar () ) f = s == '-' ? -f : f;
	for ( ; '0' <= s && s <= '9'; s = getchar () ) x = x * 10 + ( s ^ '0' );
	return x * f;
}

template<typename Tp>
inline void wint ( Tp x ) {
	if ( x < 0 ) putchar ( '-' ), x = -x;
	if ( 9 < x ) wint ( x / 10 );
	putchar ( x % 10 ^ '0' );
}

const int MAXN = 5e6, MOD = 1e9 + 7;
int n, m, K;
int pn, pr[MAXN + 5], mu[MAXN + 5], pwr[MAXN + 5], idm[MAXN + 5];
bool vis[MAXN + 5];

inline int imin ( const int a, const int b ) { return a < b ? a : b; }
inline int mul ( const long long a, const int b ) { return a * b % MOD; }
inline int sub ( int a, const int b ) { return ( a -= b ) < 0 ? a + MOD : a; }
inline int add ( int a, const int b ) { return ( a += b ) < MOD ? a : a - MOD; }
inline int mpow ( int a, int b ) {
	int ret = 1;
	for ( ; b; a = mul ( a, a ), b >>= 1 ) ret = mul ( ret, b & 1 ? a : 1 );
	return ret;
}

inline void sieve ( const int n ) {
	pwr[1] = mu[1] = idm[1] = 1;
	rep ( i, 2, n ) {
		if ( !vis[i] ) {
			mu[pr[++pn] = i] = MOD - 1;
			pwr[i] = mpow ( i, K );
			idm[i] = add ( mu[i], pwr[i] );
		}
		for ( int j = 1, t; ( t = i * pr[j] ) <= n; ++j ) {
			vis[t] = true, pwr[t] = mul ( pwr[i], pwr[pr[j]] );
			if ( !( i % pr[j] ) ) {
				idm[t] = mul ( pwr[pr[j]], idm[i] );
				break;
			}
			mu[t] = ( MOD - mu[i] ) % MOD;
			idm[t] = mul ( idm[i], idm[pr[j]] );
		}
	}
	rep ( i, 1, n ) idm[i] = add ( idm[i], idm[i - 1] );
}

int main () {
	int T = rint (); K = rint ();
	sieve ( MAXN );
	while ( T-- ) {
		n = rint (), m = rint ();
		int ans = 0;
		for ( int l = 1, r; l <= n && l <= m; l = r + 1 ) {
			r = imin ( n / ( n / l ), m / ( m / l ) );
			ans = add ( ans, mul ( mul ( n / l, m / l ),
				sub ( idm[r], idm[l - 1] ) ) );
		}
		wint ( ans ), putchar ( '\n' );
	}
	return 0;
}