P4449 于神之怒加强版 题解

原文创建时间：2021-02-28 23:53:39

我来水一篇。

题目链接

P4449

题意

给定 $k$，询问 $T$ 次，每次询问给出 $n,m$，求
$$\displaystyle \sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k$$

答案对 $10^9+7$ 取模。

解

大力推式子（认为 $n \leq m$）：

$\displaystyle \ \ \ \ \sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k$

$\displaystyle = \sum_{d=1}^nd^k\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]$

$\displaystyle = \sum_{d=1}^nd^k\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[i \perp j]$

$\displaystyle = \sum_{d=1}^nd^k\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{q|\gcd(i,j)}\mu(q)$

$\displaystyle = \sum_{d=1}^nd^k\sum_{q=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(q)\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[i|q]\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[j|q]$

$\displaystyle = \sum_{d=1}^nd^k\sum_{q=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\mu(q)\left\lfloor\frac{n}{dq}\right\rfloor\left\lfloor\frac{m}{dq}\right\rfloor$

$\displaystyle = \sum_{t=1}^n\sum_{d|t}d^k\mu\left(\frac{t}{d}\right)\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor$

$\displaystyle = \sum_{t=1}^n\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor\sum_{d|t}d^k\mu\left(\frac{t}{d}\right)$

记 $f_k(n) = n^k$，则原式

$\displaystyle = \sum_{t=1}^n\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor(\mu \ast f_k)(t)$

考虑如何筛出 $g_k = \mu \ast f_k$：

对于质数 $p$，

$$\displaylines{\begin{cases}
g_k(1) = 1 \\
g_k( p ) = p^k - 1 \\
g_k(p^\alpha) = p^{(\alpha-1)k}(p^k - 1)
\end{cases}}$$

从而线性筛解决。

质数的 $k$ 次幂，在线性筛时求出即可。

时间复杂度：$\displaystyle \mathrm{O} \left(N\log k + T\sqrt{n}\right)$

代码

#include<bits/stdc++.h>
#define REG register

using namespace std;

namespace ACAH
{
	typedef long long ll;
	
	const int bd = 5e6;
	const int maxn = 5e6 + 7;
	const int Mod = 1e9 + 7;
	
	int T, k;
	int N, M;
	int pr[maxn], pc;
	bool ip[maxn];
	ll pw[maxn], h[maxn];
	ll ans;
	
	inline ll Min(ll a, ll b) {return a < b ? a : b;}
	
	template<typename T> inline void qread(T &x)
	{
		x = 0; REG char ch = getchar();
		while(!isdigit(ch)) ch = getchar();
		while(isdigit(ch)) x = (x << 1) + (x << 3) + (ch ^ 48), ch = getchar();
	}
	
	inline ll qpow(ll a, int x)
	{
		REG ll t = a; --x;
		while(x) {if(x & 1) a = a * t % Mod; t = t * t % Mod; x >>= 1;}
		return a;
	}
	
	void linear_sieve()
	{
		h[1] = 1;
		for(REG int i = 2; i <= bd; i++) {
			if(!ip[i]) pr[++pc] = i, pw[pc] = qpow(i, k), h[i] = pw[pc] - 1;
			for(REG int j = 1; j <= pc && 1ll * i * pr[j] <= bd; j++) {
				REG int c = i * pr[j]; ip[c] = true;
				if(i % pr[j]) h[c] = h[i] * (pw[j] - 1) % Mod;
				else {h[c] = h[i] * pw[j] % Mod; break;}
			}
		}
		for(REG int i = 2; i <= bd; i++) h[i] += h[i - 1], h[i] %= Mod;
	}
	
	int work()
	{
		qread(T), qread(k);
		linear_sieve();
		while(T--) {
			qread(N), qread(M); if(N > M) swap(N, M);
			ans = 0;
			for(REG ll l = 1, r; l <= N; l = r + 1) {
				r = Min(N / (N / l), M / (M / l));
				ans += (N / l) * (M / l) % Mod * (h[r] - h[l - 1]) % Mod;
				ans %= Mod;
			}
			printf("%lld\n", (ans + Mod) % Mod);
		}
		return 0;
	}
}

int main() {return ACAH::work();}