我们先来化式子:
[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;
}