LOJ 2085: 洛谷 P1587: bzoj 4652: 「NOI2016」循环之美

题目传送门：LOJ #2085。

两个月之前做的傻题，还是有必要补一下博客。

题意简述：

求分子为不超过 (n) 的正整数，分母为不超过 (m) 的正整数的所有互不相等的分数中，有多少在 (k) 进制下的纯循环小数。

题解：

设分子为 (x)，分母为 (y)。

首先，因为要求的是互不相等的分数，取最简分数，即 (xperp y)。

其次，要求是纯循环小数，考虑竖式除法的过程，可以发现 (displaystylefrac{x}{y}) 在 (k) 进制下纯循环相当于存在正整数 (l) 使得 (xequiv xcdot k^lpmod{y})。

由于 (xperp y)，两边约去 (x) 得到 (k^lequiv 1pmod{y})，显然当 (k) 属于 (y) 的缩系中时可能成立，即 (yperp k)。

综上，答案为 (displaystylesum_{i=1}^{n}sum_{j=1}^{m}[iperp j][jperp k])。

为了方便，以下用 (adiv b) 表示 (displaystyleleftlfloorfrac{a}{b} ight floor)。答案为：

[egin{aligned}mathbf{Ans}&=sum_{i=1}^{n}sum_{j=1}^{m}[iperp j][jperp k]\&=sum_{i=1}^{n}sum_{j=1}^{m}sum_{d|i,d|j}mu(d)[jperp k]\&=sum_{d=1}^{min(n,m)}mu(d)sum_{i=1}^{ndiv d}sum_{j=1}^{mdiv d}[jdperp k]\&=sum_{d=1}^{min(n,m)}mu(d)(ndiv d)sum_{j=1}^{mdiv d}[jperp k][dperp k]\&=sum_{d=1}^{min(n,m)}mu(d)[dperp k](ndiv d)S_{[xperp k]}(mdiv d)end{aligned} ]

其中 (S_{f}(n)) 表示 (displaystylesum_{i=1}^{n}f(i))。
(S_{[xperp k]}(n)=(ndiv k)varphi(k)+S_{[xperp k]}(nmod k)) 可以 (mathcal{O}(k)) 预处理，(mathcal{O}(1)) 回答询问。

对外层 (ndiv d) 和 (mdiv d) 进行整除分块，问题转化为计算 (displaystylesum_{i=1}^{n}mu(i)[iperp k])。

设 (displaystyle S(n,k)=sum_{i=1}^{n}mu(i)[iperp k])，则有：

[egin{aligned}S(n,k)&=sum_{i=1}^{n}mu(i)[iperp k]\&=sum_{i=1}^{n}mu(i)sum_{d|i,d|k}mu(d)\&=sum_{d|k}mu(d)sum_{i=1}^{ndiv d}mu(id)\&=sum_{d|k}mu(d)sum_{i=1}^{ndiv d}mu(i)mu(d)[iperp d]\&=sum_{d|k}mu^2(d)sum_{i=1}^{ndiv d}mu(i)[iperp d]\&=sum_{d|k}mu^2(d)S(ndiv d,d)end{aligned} ]

递归，记忆化搜索即可。边界：(S(0,k)=0) 和 (displaystyle S(n,1)=sum_{i=1}^{n}mu(i)) 使用杜教筛计算。

复杂度大约为 (mathcal{O}left(n^{2/3}+sigma_0(k)sqrt{n}+k ight))。

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <map>

typedef long long LL;
const int MK = 2005;
const int S = 31622;
const int MN23 = 1000005;
const int MP = 78505;
const int MD = 25;

int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }

bool ip[MN23];
int p[MP], pc;
int mu[MN23], Smu[MN23];

inline void Init(int N) {
	mu[1] = 1;
	for (int i = 2; i <= N; ++i) {
		if (!ip[i]) p[++pc] = i, mu[i] = -1;
		for (int j = 1; j <= pc && p[j] * i <= N; ++j) {
			ip[p[j] * i] = 1;
			if (i % p[j]) mu[p[j] * i] = -mu[i];
			else break;
		}
	}
	for (int i = 1; i <= N; ++i) Smu[i] = Smu[i - 1] + mu[i];
}

int N, M, K, N23;
int A[MK], Vl[MD], cd;

std::map<int, LL> mp[MD];

LL Sum(int N, int K) {
	if (!N) return 0;
	if (K == 1 && N <= N23) return Smu[N];
	if (mp[K].count(N)) return mp[K][N];
	if (K > 1) {
		LL Ans = 0;
		for (int j = 1; j <= K; ++j)
			if (Vl[K] % Vl[j] == 0 && mu[Vl[j]])
				Ans += Sum(N / Vl[j], j);
		return mp[K][N] = Ans;
	}
	LL Ans = 1;
	for (int i = 2, j; i <= N; i = j + 1) {
		j = N / (N / i);
		Ans -= (j - i + 1) * Sum(N / i, 1);
	}
	return mp[1][N] = Ans;
}

LL Ans;

int main() {
	scanf("%d%d%d", &N, &M, &K);
	for (int i = 1; i <= K; ++i) A[i] = A[i - 1] + (gcd(i, K) == 1);
	Init(N23 = std::max((int)pow(N, 2./3), K));
	for (int i = 1; i <= K; ++i) if (K % i == 0 && (mu[i] || i == K)) Vl[++cd] = i;
	for (int i = 1, kN, kM, j; i <= N && i <= M; i = j + 1) {
		kN = N / i, kM = M / i;
		j = std::min(N / kN, M / kM);
		Ans = Ans + kN * ((LL)kM / K * A[K] + A[kM % K]) * (Sum(j, cd) - Sum(i - 1, cd));
	}
	printf("%lld
", Ans);
	return 0;
}