题面
题解
前面的题目转化可以看这篇 blog，这里直接给出式子：
[sum_{i = 1}^n sum_{j = 1}^m [i perp j][j perp k]
]
其中 ([n perp m] = [gcd(n, m) = 1])。
方法一：化简 ([i perp j])
[egin{aligned}
&sum_{i = 1}^n sum_{j = 1}^m [i perp j][j perp k]\
=&sum_{i = 1}^n sum_{j = 1}^m [j perp k] sum_{d | i, d | j} mu(d)\
=&sum_{d = 1}^n mu(d) sum_{d | i} sum_{d | j} [j perp k]\
=&sum_{d = 1}^n [d perp k] mu(d) leftlfloor frac nd
ight
floor sum_{j = 1}^{m/d} [j perp k]\
end{aligned}
]
设 (f(n) = sum_{i = 1}^n [i perp k], g(n, k) = sum_{i=1}^n mu(i)[i perp k])。
首先可以得出：
[f(n) = leftlfloor frac nk
ight
floor f(k) + f(n mod k)
]
证明可以见 M_sea 的 blog。
接下来要求出 (g(n, k))。当 (k eq 1) 时：
[egin{aligned}
g(n, k) &= sum_{i = 1}^n mu(i)[i perp k]\
&= sum_{i=1}^n mu(i) sum_{d|i, d|k} mu(d)\
&= sum_{d|k} mu(d)sum_{d|i} mu(i)\
&= sum_{d|k} mu(d)sum_{i=1}^{n/d} mu(id)\
&= sum_{d|k} mu^2(d)sum_{i=1}^{n/d} mu(i)[i perp d]\
&= sum_{d|k} mu^2(d) gleft(leftlfloor frac nd
ight
floor, d
ight)
end{aligned}
]
其中 (mu^2(x) = mu(x)^2)。
否则 (g(n, 1)) 为 (mu) 的前缀和，杜教筛即可。
方法二：化简 ([j perp k])
设 (f(n, m, k) = sum_{i=1}^nsum_{j=1}^m[i perp j][j perp k])，那么：
[egin{aligned}
f(n,m,k) &= sum_{i=1}^nsum_{j=1}^m [iperp j] sum_{d|j,d|k} mu(d)\
&= sum_{d|k}mu(d)sum_{d|j}sum_{i=1}^n[i perp j]\
&= sum_{d|k}mu(d)sum_{j=1}^{m/d}sum_{i=1}^n[iperp j][i perp d]\
&= sum_{d|k}mu(d)fleft(leftlfloor frac md
ight
floor, n, d
ight)
end{aligned}
]
边界为 (f(0, m, k) = f(n, 0, k) = 0) 以及 (f(n, m, 1) = sum_{i=1}^nsum_{j=1}^m[iperp j] = sum_{d=1}^n mu(d) lfloor frac nd floor lfloor frac md floor)。
同样可以杜教筛 (mu) 解决问题。
代码
方法一
#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#include<map>
#define RG register
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define clear(x, y) memset(x, y, sizeof(x))
inline int read()
{
int data = 0, w = 1; char ch = getchar();
while(ch != '-' && (!isdigit(ch))) ch = getchar();
if(ch == '-') w = -1, ch = getchar();
while(isdigit(ch)) data = data * 10 + (ch ^ 48), ch = getchar();
return data * w;
}
const int N(10000000), maxn(N + 10), maxk(5010);
int n, m, K, f[maxk], not_prime[maxn];
int prime[maxn], mu[maxn], smu[maxn], cnt;
std::map<std::pair<int, int>, int> sum;
void Init()
{
for(RG int i = 1; i <= K; i++) f[i] = f[i - 1] + (std::__gcd(i, K) == 1);
not_prime[1] = mu[1] = 1;
for(RG int i = 2; i <= N; i++)
{
if(!not_prime[i]) prime[++cnt] = i, mu[i] = -1;
for(RG int j = 1; j <= cnt && 1ll * i * prime[j] <= N; j++)
{
not_prime[i * prime[j]] = 1;
if(i % prime[j]) mu[i * prime[j]] = -mu[i];
else break;
}
}
for(RG int i = 1; i <= N; i++) smu[i] = smu[i - 1] + mu[i];
}
int Sf(int x) { return (x / K) * f[K] + f[x % K]; }
int Sum(int x, int k)
{
if(x == 0) return 0;
if(k == 1 && x <= N) return smu[x];
if(sum.find(std::make_pair(x, k)) != sum.end())
return sum[std::make_pair(x, k)];
int &res = sum[std::make_pair(x, k)]; res = 0;
if(k == 1)
{
res = 1; for(RG int i = 2, j; i <= x; i = j + 1)
j = x / (x / i), res -= (j - i + 1) * Sum(x / i, 1);
}
else
{
for(RG int i = 1; i * i <= k; i++) if(k % i == 0)
{
if(mu[i]) res += Sum(x / i, i);
if(i * i != k && mu[k / i]) res += Sum(x / (k / i), k / i);
}
}
return res;
}
int main()
{
#ifndef ONLINE_JUDGE
file(cpp);
#endif
n = read(), m = read(), K = read();
Init(); long long ans = 0;
for(RG int i = 1, j, min = std::min(n, m); i <= min; i = j + 1)
j = std::min(n / (n / i), m / (m / i)), ans += 1ll *
(Sum(j, K) - Sum(i - 1, K)) * (n / i) * Sf(m / i);
printf("%lld
", ans);
return 0;
}
方法二
#include <cstdio>
#include <algorithm>
#include <vector>
#include <map>
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
const int N(2010), M(1e6 + 10), LIM(1e6);
std::map<std::pair<std::pair<int, int>, int>, long long> f;
int n, m, K, mu[M], smu[M];
std::vector<int> fac[N]; std::map<int, long long> S;
long long Smu(int x)
{
if (x <= M) return smu[x];
if (S.find(x) != S.end()) return S[x];
long long &y = S[x]; y = 1;
for (int i = 2, j; i <= x; i = j + 1)
j = x / (x / i), y -= Smu(x / i) * (j - i + 1);
return y;
}
long long F(int n, int m, int K)
{
if (n == 0 || m == 0) return 0;
auto P = std::make_pair(std::make_pair(n, m), K);
if (f.find(P) != f.end()) return f[P];
long long &y = f[P]; y = 0;
if (K == 1)
{
if (n > m) std::swap(n, m);
for (int i = 1, j; i <= n; i = j + 1)
j = std::min(n / (n / i), m / (m / i)), y += (Smu(j) - Smu(i - 1)) * (n / i) * (m / i);
}
else for (auto i : fac[K]) y += F(m / i, n, i) * mu[i];
return y;
}
int main()
{
#ifndef ONLINE_JUDGE
file(cpp);
#endif
scanf("%d%d%d", &n, &m, &K), mu[1] = 1;
for (int i = 1; (i << 1) <= LIM; i++)
for (int j = (i << 1); j <= LIM; j += i) mu[j] -= mu[i];
for (int i = 1; i <= K; i++)
for (int j = i; j <= K; j += i) if (mu[i]) fac[j].push_back(i);
for (int i = 1; i <= LIM; i++) smu[i] = smu[i - 1] + mu[i];
printf("%lld
", F(n, m, K));
return 0;
}