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

link

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

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

(sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k)

(=sum_{p=1}^np^ksum_{i=1}^nsum_{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}^nmu(d)lfloorfrac n{dp} floorlfloorfrac m{dp} floor)

(=sum_{q=1}^nsum_{p|q}p^{k}mu(frac qp)lfloorfrac n{q} floorlfloorfrac m{q} floor)

注意这里求得是个数，不需要提出(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
", 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
", ans);
	}
	return 0;
}

56行，交上去一遍A