Description
给定N和M,设 p 为任一质数,求
\sum_{i=1}^{N} \sum_{j=1}^{M} [gcd(i,j) = p]
Sol
常规套路还是先把布尔表达式除一个 p
\sum_{i=1}^{N} \sum_{j=1}^{M} [gcd(\frac{i}{p},\frac{j}{p}) = 1]
然后改成枚举倍数
\sum_{i=1}^{\lfloor \frac{N}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{M}{p} \rfloor} [gcd(i,j) = 1]
[gcd(i,j)=1] 这个东西你还是很熟悉,把他莫比乌斯反演一下
\sum_{i=1}^{\lfloor \frac{N}{p} \rfloor} \sum_{j=1}^{\lfloor \frac{M}{p} \rfloor} \sum_{k|gcd(i,j)} \mu(k)
可以直接枚举质数 p
\sum_{p}^{max(N,M)} \sum_{k=1}^{\lfloor \frac{N}{p} \rfloor} \mu(k) \times \lfloor \frac{N}{kp} \rfloor \times \lfloor \frac{M}{kp} \rfloor
但是你发现在 10^7 的范围内枚举质数还 10^4 组回答不现实。再优化,考虑直接枚举 k\times p
设 T = k \times p.那么
\sum_{T=1}^{N} \lfloor \frac{N}{T} \rfloor \lfloor \frac{M}{T} \rfloor \sum_{p|T \& p \in P} \mu(\frac{T}{p})
其中 P 表示质数集
现在就做完了
Code
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int SIZE = 1e7 + 5;
int T, n, m, tot;
int isPrime[SIZE], prime[SIZE], mu[SIZE], ans[SIZE], sum[SIZE];
namespace ae86 {
const int bufl = 1 << 15;
char buf[bufl], *s = buf, *t = buf;
inline int fetch() {
if (s == t) { t = (s = buf) + fread(buf, 1, bufl, stdin); if (s == t) return EOF; }
return *s++;
}
inline int read() {
int a = 0, b = 1, c = fetch();
while (!isdigit(c))b ^= c == '-', c = fetch();
while (isdigit(c)) a = a * 10 + c - 48, c = fetch();
return b ? a : -a;
}
}
using ae86::read;
inline void getMu()
{
int siz = SIZE - 5;
isPrime[1] = mu[1] = 1;
for (register int i = 2; i <= siz; ++ i) {
if (!isPrime[i]) {
prime[++ tot] = i, mu[i] = -1;
}
for (register int j = 1; j <= tot && i * prime[j] <= siz; ++ j) {
isPrime[i * prime[j]] = 1;
if (i % prime[j]) mu[i * prime[j]] = -mu[i];
else {
mu[i * prime[j]] = 0;
break;
}
}
}
}
signed main()
{
// freopen("code.in", "r", stdin);
// freopen("code.out", "w", stdout);
T = read(); int siz = SIZE - 5;
getMu();
for (register int i = 1; i <= tot; ++ i) {
for (int j = prime[i]; j <= siz; j += prime[i]) ans[j] += mu[j / prime[i]];
}
for (register int i = 1; i <= siz; ++ i) sum[i] = sum[i - 1] + ans[i];
while (T --) {
n = read(), m = read();
int range = std::min(n, m), res = 0;
for (register int l = 1, r = 0; l <= range; l = r + 1) {
r = std::min(range, std::min(n / (n / l), m / (m / l)));
res += (n / l) * (m / l) * (sum[r] - sum[l - 1]);
}
printf("%lld\n", res);
}
return 0;
}
Comments | NOTHING