【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])

[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;
}