Description
给定a,b,d, 求
\sum_{i=1}^{a} \sum_{j=1}^{b} [gcd(i,j) = d]
Sol
还是老套路,先直接把中括号的东西除一个 d
\sum_{i=1}^{a} \sum_{j=1}^{b} [gcd(\frac{i}{d},\frac{j}{d}) = 1]
然后改成枚举约数的
\sum_{i=1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{b}{d} \rfloor} [gcd(i,j) = 1]
你可能会想把后面那个 [gcd(i,j)=1] 改成 \varphi 函数,但事实上对这个题一点帮助都没有
考虑使用莫比乌斯反演最常用的性质: \sum_{d|n} \mu(d) = [n=1]
相当于上面这个式子的 n 就是 gcd(i,j) ,枚举 k | gcd(i,j) ,对 \mu(k) 求和
\sum_{i=1}^{\lfloor \frac{a}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{b}{d} \rfloor} \sum_{k|gcd(i,j)} \mu(k)
然后改成枚举 gcd(i,j) 的约数
\sum_{k=1}^{\lfloor \frac{a}{d} \rfloor} \sum_{i=1}^{\lfloor {\frac{a}{dk} \rfloor}} \sum_{j=1}^{\lfloor \frac{b}{dk} \rfloor} \mu(k)
事实上第一个和式 \sum_{k=1}^{\lfloor \frac{a}{d} \rfloor} 的上届应该是 max(\lfloor \frac{a}{d} \rfloor, \lfloor \frac{b}{d}\rfloor),但是那个上面的位置太小了,写在这里
后面那一对东西可以直接整出分块算掉,所以最终的式子是:
\sum_{k=1}^{\lfloor \frac{a}{d} \rfloor} \mu(k) \times \lfloor {\frac{a}{dk} \rfloor} \times \lfloor \frac{b}{dk} \rfloor
Code
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int SIZE = 5e4 + 5;
int T, a, b, d, tot;
int prime[SIZE], isPrime[SIZE], mu[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 = 5e4;
isPrime[1] = mu[1] = 1;
for (int i = 2; i <= siz; ++ i) {
if (!isPrime[i]) {
prime[++ tot] = i, mu[i] = -1;
}
for (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();
getMu();
for (int i = 1; i <= 5e4; ++ i) sum[i] = sum[i - 1] + mu[i];
while (T --) {
a = read(), b = read(), d = read();
a /= d, b /= d;
int range = std::min(a, b), ans = 0;
for (int l = 1, r; l <= range; l = r + 1) {
if (!(range / l)) r = range;
else r = std::min(range, std::min(a / (a / l), b / (b / l)));
ans += (sum[r] - sum[l - 1]) * (a / l) * (b / l);
}
printf("%lld\n", ans);
}
return 0;
}
Comments | NOTHING