【题解】Luogu-P4449 于神之怒加强版

P4449 于神之怒加强版

Description

• 多测，数据组数为 $T$。

• 给定 常数 $k$ 和整数 $n,m$，计算

$$\left[\sum _{i=1}^n\sum _{j=1}^m\gcd (i,j{)}^k\right]\mathrm{m}\mathrm{o}\mathrm{d}(10^9+7)$$

• 对于全部的测试点，保证 $1\le T\le 2\times 10^3,1\le n,m,k\le 5\times 10^6$。

Solution

不妨设 $n\le m$。
$$\begin{array}{rl}\sum _{i=1}^n\sum _{j=1}^m\gcd (i,j{)}^k& =\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m{d}^k[\gcd (i,j)=d]\\ & =\sum _{d=1}^n{d}^k\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[\gcd (i,j)=1]\\ & =\sum _{d=1}^n{d}^k\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\sum _{p\mid \gcd (i,j)}\mu (p)\\ & =\sum _{d=1}^n{d}^k\sum _{p=1}^n\mu (p)\lfloor \frac{n}{dp}\rfloor \lfloor \frac{m}{dp}\rfloor \\ & =\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor \sum _{p\mid T}\mu (p){\left(\frac{T}{p}\right)}^k\\ & =\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor (\mu \ast {\mathrm{Id}}_k)(T)\end{array}$$
设 $f(n)=(\mu \ast {\mathrm{Id}}_k)(n)$，由 $\mu$ 与 ${\mathrm{Id}}_k$ 都是积性函数可得 $f$ 是积性函数。

考虑 $f$ 的线性筛。

• $i$ 为质数：$f(i)=i^k-1$
• $p_j\nmid i$：$f(i\cdot p_j)=f(i)f(p_j)$
• $p_j\mid i$：$f(i\cdot p_j)=f(i)p_j^k$

预处理 $\mathrm{O}(n)$，整除分块 $\mathrm{O}(\sqrt{n})$，总时间复杂度为 $\mathrm{O}(n+T\sqrt{n})$。

Code

// 18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

const int MAXN = 5e6 + 5;
const int N = 5e6;
const int MOD = 1e9 + 7;

int qpow(int a, int b)
{
	int base = a, ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = (ll)ans * base % MOD;
		}
		base = (ll)base * base % MOD;
		b >>= 1;
	}
	return ans;
}

int p[MAXN], xsy[MAXN], f[MAXN], sum[MAXN];
bool vis[MAXN];

void pre(int k)
{
	f[1] = sum[1] = 1;
	for (int i = 2; i <= N; i++)
	{
		if (!vis[i])
		{
			p[++p[0]] = i;
			xsy[i] = qpow(i, k);
			f[i] = (xsy[i] - 1 + MOD) % MOD;
		}
		for (int j = 1; j <= p[0] && i * p[j] <= N; j++)
		{
			vis[i * p[j]] = true;
			if (i % p[j] == 0)
			{
				f[i * p[j]] = (ll)f[i] * xsy[p[j]] % MOD;
				break;
			}
			f[i * p[j]] = (ll)f[i] * f[p[j]] % MOD;
		}
		sum[i] = (sum[i - 1] + f[i]) % MOD;
	}
}

int getsum(int l, int r)
{
	return (sum[r] - sum[l - 1] + MOD) % MOD;
}

int block(int n, int m)
{
	int res = 0;
	for (int l = 1, r; l <= n; l = r + 1)
	{
		int k1 = n / l, k2 = m / l;
		r = min(n / k1, m / k2);
		res = (res + (ll)k1 * k2 % MOD * getsum(l, r) % MOD) % MOD;
	}
	return res;
}

int main()
{
	int t, k;
	scanf("%d%d", &t, &k);
	pre(k);
	while (t--)
	{
		int n, m;
		scanf("%d%d", &n, &m);
		if (n > m)
		{
			swap(n, m);
		}
		printf("%d\n", block(n, m));
	}
	return 0;
}