洛谷 P4449 于神之怒加强版

题意

求

[sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j)^k(mod 1e9+7) ]

思路

还是直接淦式子

[egin{align*}&sum_{i=1}^{n}sum_{j=1}^{m}gcd(i,j)^k\=&sum_{i=1}^{n}sum_{j=1}^{m}sum_{d=1}^{min(n,m)}d^k[gcd(i,j)=d]\=&sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{lfloorfrac nd floor}sum_{j=1}^{lfloorfrac md floor}[gcd(i,j)=1]\=&sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{lfloorfrac nd floor}sum_{j=1}^{lfloorfrac md floor}sum_{x|gcd(i,j)}mu(x)\=&sum_{d=1}^{min(n,m)}d^ksum_{i=1}^{lfloorfrac nd floor}sum_{j=1}^{lfloorfrac md floor}sum_{x|i,x|j}mu(x)\=&sum_{d=1}^{min(n,m)}d^ksum_{x=1}^{min(lfloorfrac nd floor,lfloorfrac md floor)}mu(x)sum_{i=1}^{lfloorfrac nd floor}[x|i]sum_{j=1}^{lfloorfrac md floor}[x|j]\=&sum_{d=1}^{min(n,m)}d^ksum_{x=1}^{min(lfloorfrac nd floor,lfloorfrac md floor)}mu(x)lfloorfrac n{dx} floorlfloorfrac m{dx} floorend{align*} ]

令(P=dx)，则原式等于

[sum_{P=1}^{min(n,m)}lfloorfrac n{P} floorlfloorfrac m{P} floorsum_{d|P}d^kmu(frac Pd) ]

显然前面的(lfloorfrac n{P} floorlfloorfrac m{P} floor)部分可以分块求解。

现在考虑后面的一部分，令

[g(n)=sum_{d|n}d^kmu(frac nd) ]

容易得出这个函数是积性函数，所以我们就可以线性筛然后求出其前缀和

然后就做完了

代码

/*
Author:loceaner
莫比乌斯反演
*/
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define int long long
using namespace std;

const int A = 5e6 + 11;
const int B = 1e6 + 11;
const int mod = 1e9 + 7;
const int inf = 0x3f3f3f3f;

inline int read() {
	char c = getchar();
	int x = 0, f = 1;
	for( ; !isdigit(c); c = getchar()) if(c == '-') f = -1;
	for( ; isdigit(c); c = getchar()) x = x * 10 + (c ^ 48);
	return x * f;
}

bool vis[A];
int T, n, m, k, f[A], g[A], p[A], cnt, sum[A];

inline int power(int a, int b) {
	int res = 1;
	while (b) {
		if (b & 1) res = res * a % mod;
		a = a * a % mod, b >>= 1;
	}
	return res;
}

inline int mo(int x) {
	if(x > mod) x -= mod;
	return x;
}

inline void work() {
	g[1] = 1;
	int maxn = 5e6 + 1;
	for (int i = 2; i <= maxn; i++) {
		if (!vis[i]) { p[++cnt] = i, f[cnt] = power(i, k), g[i] = mo(f[cnt] - 1 + mod); }
		for (int j = 1; j <= cnt && i * p[j] <= maxn; j++) {
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) { g[i * p[j]] = g[i] * 1ll * f[j] % mod; break; }
			g[i * p[j]] = g[i] * 1ll * g[p[j]] % mod;
		}
	}
	for (int i = 2; i <= maxn; i++) g[i] = (g[i - 1] + g[i]) % mod;
}

inline int abss(int x) {
	while (x < 0) x += mod;
	return x;
}

signed main() {
	T = read(), k = read();
	work();
	while (T--) {
		n = read(), m = read();
		int maxn = min(n, m), ans = 0;
		for (int l = 1, r; l <= maxn; l = r + 1) {
			r = min(n / (n / l), m / (m / l));
			(ans += abss(g[r] - g[l - 1]) * 1ll * (n / l) % mod * (m / l) % mod) %= mod;
		}
		ans = (ans % mod + mod) % mod;
		cout << ans << '
';
	}
	return 0;
}