[BZOJ4407] 于神之怒加强版

---

Portal

我们先来化式子:

[Ans = sum_{i = 1}^{n} sum_{j = 1}^{m} (i, j) ^{k} pmod {1e9 + 7}\ ]

首先是套路枚举GCD + 反演

[Ans = sum_{i}sum_{j}sum_{l = 1} ^ {min(n, m)} l ^ {k} ~cdot~ [(i, j) == l] \ = sum_{l = 1} ^ {min(n, m)} l ^ {k} sum_{i}sum_{j} [(frac{i}{l}, frac{j}{l})= 1] \ = sum_{l = 1} ^ {min(n, m)} l ^ {k} sum_{i}^{frac{n}{l}}sum_{j}^{frac{m}{l}} [(i,j)= 1] \ = sum_{l = 1} ^ {min(n, m)} l ^ {k} sum_{i}sum_{j} sum_{d} [d | i][d | j] mu(d) \ = sum_{l}l^k sum_{d} mu(d)frac{n}{ld}frac{m}{ld} \ ]

然后就是合并(sum)的套路: 观察到式子中有两个sum结合起来的部分.我们可以枚举结合部,然后再枚举另一部.

[Ans = sum_{l}sum_{d} l^kmu(d)frac{n}{ld}frac{m}{ld}\ =sum_{c}sum_{p | c} p ^ k frac{}{} mu(frac{c}{p}) frac{}{} frac{n}{c} frac{m}{c} ]

定义$$f(n) = sum_{d | n} d ^ k frac{}{} mu(frac{n}{d}) $$,发现这就是(Id * mu),那么f为积性函数, 可以采用Min25线性筛筛一下前缀和.然后整除分块就做完了.

---

以下是线筛筛积性函数的通用办法.

明显的, 每一个数都只会被它自己的最小素因子筛掉, 那么我们记录一个mp[i]表示i质因数分解中的(p_1 ^ {a_1})

那么就很好做了:

对于某一个(prime_j):

((1)) 如果 i % prime[j] != 0, F[i * prime[j]] = F[i] * F[prime[j]]

否则((2)) 当mp[i] == i, 表示i是一个质数的幂,我们可以通过递推或其他操作得到.否则 F[i * prime[j]] = F[i / mp[i]] * F[mp[i] * prime[j]]

这样的话,把最小质因子除掉了,他们就是互质的.

Code

#include<bits/stdc++.h>
using namespace std;
#define rep(i, a, b) for(int i = (a), i##_end_ = (b); i <= i##_end_; ++i)
#define drep(i, a, b) for(int i = (a), i##_end_ = (b); i >= i##_end_; --i)
#define clar(a, b) memset((a), (b), sizeof(a))
#define debug(...) fprintf(stderr, __VA_ARGS__)
typedef long long LL;
typedef long double LD;
int read() {
    char ch = getchar();
    int x = 0, flag = 1;
    for (;!isdigit(ch); ch = getchar()) if (ch == '-') flag *= -1;
    for (;isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
    return x * flag;
}
void write(int x) {
    if (x < 0) putchar('-'), x = -x;
    if (x >= 10) write(x / 10);
    putchar(x % 10 + 48);
}

const int Mod = 1e9 + 7;
const int Maxn = 5e6 + 9;
int fpm(int base, int tims) {
	int r = 1; base %= Mod;
	while (tims) {
		if (tims & 1) r = 1ll * base * r % Mod;
		base = 1ll * base * base % Mod;
		tims >>= 1;
	}
	return r;
}

int T, k;

int prime[Maxn], isnprime[Maxn], mu[Maxn], tot, F[Maxn], low[Maxn], pk[Maxn];
int sumfixF[Maxn];
inline int getPk(int u) { return pk[u] ? pk[u] : (pk[u] = fpm(u, k)); }
void linearSieve() {
	mu[1] = 1; F[1] = 1;
	rep (i, 2, Maxn - 1) {
		if (!isnprime[i]) prime[++tot] = i, mu[i] = -1, F[i] = (getPk(i) - 1ll + Mod) % Mod, low[i] = i;
		for (int k, j = 1; j <= tot && (k = prime[j] * i) < Maxn; ++j) {
			isnprime[k] = 1;
			if (i % prime[j] == 0) {
				low[k] = low[i] * prime[j], mu[k] = 0;
				if (low[k] == k) F[k] = 1ll * F[i] * getPk(prime[j]) % Mod; 
/**/			else F[k] = 1ll * F[k / low[k]] * F[low[k]] % Mod;
				break;
			} else {
				mu[k] = mu[i] * mu[prime[j]];
				F[k] = 1ll * F[i] * F[prime[j]] % Mod;
				low[k] = prime[j];
			}
		}
	}

	rep (i, 1, Maxn - 1) sumfixF[i] = (1ll * sumfixF[i - 1] + F[i]) % Mod; 
	
}

vector <pair<int, int> > s;
void init() {
	T = read(), k = read();
	rep (i, 1, T) {
		int u = read(), v = read();
		s.push_back(make_pair(u, v));
	}
	
	linearSieve();
}

void solve() {
	rep (i, 1, T) {
		int n = s[i - 1].first, m = s[i - 1].second;
		int Limit = min(n, m); LL ans = 0;
		for (int r, l = 1; l <= Limit; l = r + 1) {
			r = min(min(n / (n / l), m / (m / l)), Limit);
			(ans += 1ll * (1ll * sumfixF[r] - sumfixF[l - 1] + Mod) % Mod * (n / l) % Mod * (m / l) % Mod) %= Mod;
		}
		printf("%d
", ans);
	}
}

int main() {
	freopen("BZOJ4407.in", "r", stdin);
	freopen("BZOJ4407.out", "w", stdout);

	init();
	solve();

#ifdef Qrsikno
    debug("
Running time: %.3lf(s)
", clock() * 1.0 / CLOCKS_PER_SEC);
#endif
    return 0;
}