P4449 於神之怒加強版

Description

• 多測，資料組數為 \(T\)。

• 給定 常數 \(k\) 和整數 \(n, m\)，計算

\[\left[\sum_{i = 1}^n \sum_{j = 1}^m \gcd(i, j)^k \right] \bmod (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{aligned} \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}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum_{j = 1}^{\left\lfloor\frac{m}{d}\right\rfloor} [\gcd(i, j) = 1] \\ & = \sum_{d = 1}^n d^k \sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor} \sum_{j = 1}^{\left\lfloor\frac{m}{d}\right\rfloor} \sum_{p\mid \gcd(i, j)} \mu(p) \\ & = \sum_{d = 1}^n d^k \sum_{p = 1}^n \mu(p) \left\lfloor\dfrac{n}{dp}\right\rfloor \left\lfloor\dfrac{m}{dp}\right\rfloor \\ & = \sum_{T = 1}^n \left\lfloor\dfrac{n}{T}\right\rfloor \left\lfloor\dfrac{m}{T}\right\rfloor \sum_{p\mid T} \mu(p) \left(\dfrac{T}{p}\right)^k \\ & = \sum_{T = 1}^n \left\lfloor\dfrac{n}{T}\right\rfloor \left\lfloor\dfrac{m}{T}\right\rfloor (\mu * \operatorname{Id}_k)(T) \end{aligned} \]

設 \(f(n) = (\mu * \operatorname{Id}_k)(n)\)

，由 \(\mu\) 與 \(\operatorname{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\)

預處理 \(\Omicron(n)\)，整除分塊 \(\Omicron(\sqrt{n})\)，總時間複雜度為 \(\Omicron(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;
}