【NOI2016】循环之美
题面
题解
前面的题目转化可以看这篇 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]\)
\[\begin{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) \left\lfloor \frac nd \right\rfloor \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) = \left\lfloor \frac nk \right\rfloor f(k) + f(n \bmod k)
\]
证明可以见 M_sea 的 blog。
接下来要求出 \(g(n, k)\)。当 \(k \neq 1\) 时：
\[\begin{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) g\left(\left\lfloor \frac nd \right\rfloor, d\right)
\end{aligned}
\]
其中 \(\mu^2(x) = \mu(x)^2\)。
否则 \(g(n, 1)\) 为 \(\mu\) 的前缀和，杜教筛即可。
upd: 突然发现自己 naive 了。
事实上 \(S(n) = \sum_{i=1}^n \mu(i) [i \perp k]\) 是积性函数的和，构造函数 \(h(n) = [n \perp k]\) 和 \(\mu(n) [n \perp k]\) 卷起来就可以杜教筛了。
\[S(n) = 1 - \sum_{i=2}^n [i \perp k] S\left(\left \lfloor \frac ni \right\rfloor\right)
\]
方法二：化简 \([j \perp k]\)
设 \(f(n, m, k) = \sum_{i=1}^n\sum_{j=1}^m[i \perp j][j \perp k]\)，那么：
\[\begin{aligned}
f(n,m,k) &= \sum_{i=1}^n\sum_{j=1}^m [i\perp j] \sum_{d|j,d|k} \mu(d)\\
&= \sum_{d|k}\mu(d)\sum_{j|d}\sum_{i=1}^n[i \perp j]\\
&= \sum_{d|k}\mu(d)\sum_{j=1}^{m/d}\sum_{i=1}^n[i\perp j][i \perp d]\\
&= \sum_{d|k}\mu(d)f\left(\left\lfloor \frac md \right\rfloor, n, d\right)
\end{aligned}
\]
边界为 \(f(0, m, k) = f(n, 0, k) = 1\) 以及 \(f(n, m, 1) = \sum_{i=1}^n\sum_{j=1}^m[i\perp j] = \sum_{d=1}^n \mu(d) \lfloor \frac nd\rfloor \lfloor \frac md\rfloor\)。
同样可以杜教筛 \(\mu\) 解决问题。
代码
方法一
这里放的是杜教筛的做法，递归做法网上应该到处都是。
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <map>
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
inline int read()
{
int data = 0, w = 1; char ch = getchar();
while (ch != '-' && (ch < '0' || ch > '9')) ch = getchar();
if (ch == '-') w = -1, ch = getchar();
while (ch >= '0' && ch <= '9') data = data * 10 + (ch ^ 48), ch = getchar();
return data * w;
}
const int N(1e6 + 10), M(2010);
int n, m, B, L, K, f[M], h[N], not_prime[N], prime[N >> 3], cnt, g[N];
inline int Sf(int x) { int y = x / K; return y * f[K] + f[x - y * K]; }
std::map<int, int> sg;
int S(int n)
{
if (n <= L) return g[n];
if (sg.find(n) != sg.end()) return sg[n];
int res = 1, _n = std::sqrt(n);
for (int i = 2; i <= _n; i++) res -= h[i] * S(n / i);
for (int i = _n - (_n * _n == n), k, t, p = Sf(_n); i >= 1; i--)
k = n / i, res -= S(i) * ((t = Sf(k)) - p), p = t;
return sg[n] = res;
}
int main()
{
#ifndef ONLINE_JUDGE
file(cpp);
#endif
n = read(), m = read(), K = read(), not_prime[1] = g[1] = h[1] = 1;
L = std::pow(B = std::min(n, m), 2. / 3);
for (int i = 1; i <= K; i++) f[i] = f[i - 1] + (std::__gcd(i, K) == 1);
for (int i = 2; i <= L; i++)
{
if (!not_prime[i]) prime[++cnt] = i, g[i] = -(h[i] = K % i != 0);
for (int j = 1; j <= cnt && i * prime[j] <= L; j++)
{
not_prime[i * prime[j]] = 1, h[i * prime[j]] = h[i] & h[prime[j]];
if (i % prime[j]) g[i * prime[j]] = g[i] * g[prime[j]];
else break;
}
}
for (int i = 2; i <= L; i++) g[i] += g[i - 1];
long long ans = 0;
for (int i = 1, j, k, l, p = 0, t; i <= B; p = t, i = j + 1)
j = std::min(n / (k = n / i), m / (l = m / i)), ans += 1ll * ((t = S(j)) - p) * k * Sf(l);
printf("%lld\n", 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\n", F(n, m, K));
return 0;
}
