luogu4449 于神之怒加强版(莫比乌斯反演)

link

给定n,m,k,计算\(\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k\)对1000000007取模的结果

多组数据，T<=2000,1<=N,M,K<=5000000
推式子

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

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

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

\(=\sum_{p=1}^np^{k}\sum_{i=1}^{n/p}\sum_{j=1}^{m/p}\sum_{d|i,d|j}\mu(d)\)

\(=\sum_{p=1}^np^{k}\sum_{d=1}^n\mu(d)\lfloor\frac n{dp}\rfloor\lfloor\frac m{dp}\rfloor\)

\(=\sum_{q=1}^n\sum_{p|q}p^{k}\mu(\frac qp)\lfloor\frac n{q}\rfloor\lfloor\frac m{q}\rfloor\)

注意这里求得是个数，不需要提出\(p^2\)和\(d^2\)，我式子推错了两次。。。

还是枚举倍数对于所有q处理\(\sum_{p|q}p^{k}\mu(\frac qp)\)，然后打数论分块

注意这里如果定义p为1e9+7就不要再用p了。。。

#include <cstdio>
#include <functional>
using namespace std;

int n, prime[5000010], mu[5000010], tot, fuck = 5000000, p = 1000000007;
int s[5000010];
bool vis[5000010];

int qpow(int x, int y)
{
	int res = 1;
	while (y > 0)
	{
		if (y & 1) res = res * (long long)x % p;
		x = x * (long long)x % p;
		y >>= 1;
	}
	return res;
}

int main()
{
	int t, k; scanf("%d%d", &t, &k);
	mu[1] = 1;
	for (int i = 2; i <= fuck; i++)
	{
		if (vis[i] == 0) prime[++tot] = i, mu[i] = -1;
		for (int j = 1; j <= tot && i * prime[j] <= fuck; j++)
		{
			vis[i * prime[j]] = true;
			if (i % prime[j] == 0) break;
			mu[i * prime[j]] = -mu[i];
		}
	}
	for (int pp = 1; pp <= fuck; pp++)
	{
		int sb = qpow(pp, k);
		for (int q = pp, d = 1; q <= fuck; q += pp, d++)
			s[q] = (s[q] + sb * mu[d]) % p;
	}
	for (int i = 1; i <= fuck; i++)
	{
		// printf("s[%d] = %d\n", i, s[i]);
		s[i] = (s[i] + s[i - 1]) % p;
	}
	while (t --> 0)
	{
		int n, m, ans = 0;
		scanf("%d%d", &n, &m); if (n > m) swap(n, m);
		for (int i = 1, j; i <= n; i = j + 1)
			j = min(n / (n / i), m / (m / i)), ans = (ans + (s[j] - s[i - 1]) * (long long)(n / i) % p * (m / i) % p) % p;
		if (ans < 0) ans += p;
		printf("%d\n", ans);
	}
	return 0;
}

56行，交上去一遍A