【BZOJ4652】【UOJ221】【NOI2016】循环之美

【题目链接】

• BZOJ
  
• UOJ
  

【思路要点】

• 通过在十进制下找规律，我们发现分数\(\frac{x}{y}\)在\(k\)进制下为纯循环小数当且仅当\(gcd(y,k)=1\)。
• 稍加分析，我们发现上面这一点并不难证明。
• 那么，原题要求的式子应当是\(\sum_{i=1}^{N}\sum_{j=1}^{M}\epsilon(gcd(i,j))*\epsilon(gcd(j,k))\)。
• \(=\sum_{i=1}^{N}\sum_{j=1}^{M}\epsilon(gcd(j,k))\sum_{d/i,d/j}\mu(d)\)
• \(=\sum_{d\bot k}\mu(d)\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}\epsilon(gcd(j,k))\)
• \(=\sum_{d\bot k}\mu(d)\lfloor\frac{N}{d}\rfloor\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}\epsilon(gcd(j,k))\)
  
• 通过适当地预处理，\(\lfloor\frac{N}{d}\rfloor\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}\epsilon(gcd(j,k))\)可以在\(O(1)\)的时间内计算得出。
• 因此，现在直接求解这个式子的时空复杂度是\(O(N)\)的，可以得到84分，是十分可观的分数。
• 现在问题在于快速如何求得\(\sum_{d\bot k}\mu(d)\)的前缀和。
• 令\(S(N)=\sum_{d\bot k}^{N}\mu(d)\)，\(M(N)=\sum_{i=1}^{N}\mu(i)\)。
• 那么有\(S(N)=M(N)-\sum_{d=2,d/k}^{N}\sum_{i=1}^{N}[gcd(i,k)=d]*\mu(i)\)。
• \(=M(N)-\sum_{d=2,d/k}^{N}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\epsilon(gcd(i,\frac{k}{d}))*\mu(id)\)
  
• \(=M(N)-\sum_{d=2,d/k}^{N}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\epsilon(gcd(i,\frac{k}{d}))*\mu(i)*\mu(d)*\epsilon(gcd(i,d))\)
  
• \(=M(N)-\sum_{d=2,d/k}^{N}\mu(d)\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\epsilon(gcd(i,k))*\mu(i)\)
  
• \(=M(N)-\sum_{d=2,d/k}^{N}\mu(d)*S(\lfloor\frac{N}{d}\rfloor)\)
  
• 那么，用杜教筛即可解决这个函数的前缀和。
• 时间复杂度\(O(\sqrt{NK}+N^{\frac{2}{3}})\)。

【代码】

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 2e6 + 5;
const int MAXK = 2005;
template <typename T> void chkmax(T &x, T y) {x = max(x, y); }
template <typename T> void chkmin(T &x, T y) {x = min(x, y); } 
template <typename T> void read(T &x) {
	x = 0; int f = 1;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
	for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
	x *= f;
}
template <typename T> void write(T x) {
	if (x < 0) x = -x, putchar('-');
	if (x > 9) write(x / 10);
	putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {
	write(x);
	puts("");
}
int n, m, k;
int pre[MAXK], pri[MAXK];
int tot, prime[MAXN], f[MAXN];
int miu[MAXN], smiu[MAXN], sum[MAXN];
int smiun[MAXN], smium[MAXN];
int sumn[MAXN], summ[MAXN];
void init() {
	miu[1] = sum[1] = smiu[1] = 1;
	for (int i = 2; i < MAXN; i++) {
		if (f[i] == 0) {
			prime[++tot] = f[i] = i;
			miu[i] = -1;
		}
		smiu[i] = smiu[i - 1] + miu[i];
		sum[i] = sum[i - 1] + miu[i] * pri[i % k];
		for (int j = 1; j <= tot && prime[j] <= f[i]; j++) {
			int tmp = i * prime[j];
			if (tmp >= MAXN) break;
			f[tmp] = prime[j];
			if (prime[j] == f[i]) miu[tmp] = 0;
			else miu[tmp] = -miu[i];
		}
	}
}
int gcd(int x, int y) {
	if (y == 0) return x;
	else return gcd(y, x % y);
}
long long func(int n, int m) {
	long long ans = m / k * pre[k] + pre[m % k];
	return ans * n;
}
int presumn(int x) {
	if (x < MAXN) return smiu[x];
	int y = n / x, ans = 1;
	if (smiun[y] != -1) return smiun[y];
	int now = 2;
	while (now <= x) {
		int nxt = x / (x / now) + 1;
		ans -= (nxt - now) * presumn(x / now);
		now = nxt;
	}
	return smiun[y] = ans;
}
int presumm(int x) {
	if (x < MAXN) return smiu[x];
	int y = m / x, ans = 1;
	if (smium[y] != -1) return smium[y];
	int now = 2;
	while (now <= x) {
		int nxt = x / (x / now) + 1;
		ans -= (nxt - now) * presumm(x / now);
		now = nxt;
	}
	return smium[y] = ans;
}
int sieven(int x) {
	if (x < MAXN) return sum[x];
	int y = n / x, ans = presumn(x);
	if (sumn[y] != -1) return sumn[y];
	for (int i = 1; i * i <= k; i++) {
		if (k % i) continue;
		int j = k / i;
		if (i >= 2 && i <= x) ans -= miu[i] * sieven(x / i);
		if (j != i && j >= 2 && j <= x) ans -= miu[j] * sieven(x / j);
	}
	return sumn[y] = ans;
}
int sievem(int x) {
	if (x < MAXN) return sum[x];
	int y = m / x, ans = presumm(x);
	if (summ[y] != -1) return summ[y];
	for (int i = 1; i * i <= k; i++) {
		if (k % i) continue;
		int j = k / i;
		if (i >= 2 && i <= x) ans -= miu[i] * sievem(x / i);
		if (j != i && j >= 2 && j <= x) ans -= miu[j] * sievem(x / j);
	}
	return summ[y] = ans;
}
int main() {
	read(n), read(m), read(k);
	for (int i = 1; i <= k; i++) {
		pri[i] = gcd(i, k) == 1;
		pre[i] = pre[i - 1] + pri[i];
	}
	init();
	int now = 1, last = 0;
	long long ans = 0;
	memset(sumn, -1, sizeof(sumn));
	memset(summ, -1, sizeof(summ));
	memset(smiun, -1, sizeof(smiun));
	memset(smium, -1, sizeof(smium));
	while (now <= n && now <= m) {
		int nxt = n / (n / now);
		int mxt = m / (m / now);
		if (nxt <= mxt) {
			ans += (sieven(nxt) - last) * func(n / now, m / now);
			now = nxt + 1; last = sieven(nxt);
		} else {
			ans += (sievem(mxt) - last) * func(n / now, m / now);
			now = mxt + 1; last = sievem(mxt);
		}
	}
	writeln(ans);
	return 0;
}