【loj

目录

• description
• solution
• accepted code
• details

---

description

牛牛是一个热爱算法设计的高中生。在他设计的算法中，常常会使用带小数的数进行计算。牛牛认为，如果在 (k) 进制下，一个数的小数部分是纯循环的，那么它就是美的。

现在，牛牛想知道：对于已知的十进制数 (n) 和 (m)，在 (k) 进制下，有多少个数值上互不相等的纯循环小数，可以用分数 (frac{x}{y}) 表示，其中 (1leq x leq n, 1 leq y leq m)，且 (x, y) 是整数。

一个数是纯循环的，当且仅当其可以写成以下形式：

[a.dot{c_1}c_2c_3dots c_{p-1}dot{c_{p}} ]

其中，(a) 是一个整数，(pgeq 1)；对于 (1leq i leq p)，(c_i) 是 (k) 进制下的一位数字。

例如，在十进制下，(0.45454545cdots=0.dot{4}dot{5}) 是纯循环的，它可以用 (frac{5}{11})、(frac{10}{22}) 等分数表示；在十进制下，(0.1666666cdots=0.1dot{6}) 则不是纯循环的，它可以用 (frac{1}{6}) 等分数表示。

需要特别注意的是，我们认为一个整数是纯循环的，因为它的小数部分可以表示成 (0) 的循环或是 (k-1) 的循环；而一个小数部分非 (0) 的有限小数不是纯循环的。

solution

考虑既约分数 (frac{j}{i}) 是纯循环的等价条件：存在 (p) 使得 (j imes k^p=j mod i)。

由于 (j, i) 互质，那么上述条件等价于 (k, i) 也互质。因此答案为：

[sum_{i=1}^{m}sum_{j=1}^{n}[gcd(i,k) = 1][gcd(i,j)=1] ]

首先套路反演：

[sum_{d|k}mu(d)sum_{p}mu(p) imeslfloorfrac{n}{p} floor imeslfloorfrac{m}{{ m lcm}(d,p)} floor ]

二度反演：

[sum_{a|b|d|k}mu(d)sum_{p'}mu(b imes p') imeslfloorfrac{n}{b imes p'} floor imeslfloorfrac{m imes a}{d imes b imes p'} floor ]

枚举四元组 ((a, b, d, k))，记 (n' = lfloorfrac{n}{b} floor, m'=lfloorfrac{m imes a}{d imes b} floor)，则我们相当于求：

[sum_{p'}mu(b imes p') imeslfloorfrac{n'}{p'} floor imeslfloorfrac{m'}{p'} floor ]

后面部分很明显可以分块。考虑怎么快速求 (S_{b}(n) = sum_{i=1}^{n}mu(b imes i))，三度反演：

[egin{aligned} sum_{i=1}^{n}mu(b imes i) &=sum_{i=1}^{n}mu(b) imesmu(i) imes[gcd(b,i)=1] \ &=sum_{d|b}mu(b) imesmu(d) imes S_{d}(lfloorfrac{n}{d} floor) end{aligned} ]

当 (b = 1) 时直接杜教筛，否则就做如上的过程。不难发现 (b ot=1) 时转移无环。

对于最后一步 (b ot = 1) 的时间复杂度分析：一共会有 (O(d(k) imessqrt{n})) 个状态，每个状态有 (O(d(k))) 种转移。所以可以通过。

accepted code

#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>
using namespace std;

typedef long long ll;

const int MAXK = 2000;
const int MAXN = 1000000;
const ll INF = ll(1E12);

int mu[MAXN + 5], prm[MAXN + 5], pcnt;
bool nprm[MAXN + 5];

vector<int>dv[MAXK + 5];
void init() {
	for(int i=1;i<=MAXK;i++)
		for(int j=i;j<=MAXK;j+=i)
			dv[j].push_back(i);
			
	mu[1] = 1;
	for(int i=2;i<=MAXN;i++) {
		if( !nprm[i] ) prm[++pcnt] = i, mu[i] = -1;
		for(int j=1;i*prm[j]<=MAXN;j++) {
			nprm[i*prm[j]] = true;
			if( i % prm[j] == 0 ) {
				mu[i*prm[j]] = 0;
				break;
			}
			else mu[i*prm[j]] = -mu[i];
		}
	}
}

ll f[20][MAXN + 5]; int id1[MAXN + 5], id2[MAXN + 5], tot;
int id(int x) {return (x <= MAXN ? id1[x] : id2[INF / x]);}

ll smu[20][MAXN + 5]; int A[MAXK + 5], B[20], cnt;
ll sum(int x, int n) {
	ll &ans = f[x][id(n)];
	if( ans != INF ) return ans;
	if( n <= MAXN ) return ans = smu[x][n];
	
	if( x == 1 ) {
		ans = 1;
		for(int i=2,j;i<=n;i=j+1) {
			int p = (n / i); j = n / p;
			ans -= (j - i + 1)*sum(x, p);
		}
	} else {
		ans = 0;
		for(unsigned i=0;i<dv[B[x]].size();i++) {
			int d = dv[B[x]][i];
			ans += mu[d]*sum(A[d], n / d);
		}
		ans *= mu[B[x]];
	}
	return ans;
}

void get_id(int n, int m) {
	tot = 0;
	for(int i=1;i<=n&&i<=m;i++) {
		i = min(n / (n / i), m / (m / i));
		if( i <= MAXN ) id1[i] = (++tot);
		else id2[INF / i] = (++tot);
		
		for(int j=1;j<=cnt;j++) f[j][tot] = INF;
	}
}
ll solve(int b, int n, int m) {
	get_id(n, m); ll ret = 0;
	for(int i=1,j;i<=n&&i<=m;i=j+1) {
		int p = (n / i), q = (m / i);
		j = min(n / p, m / q);
		ret += 1LL*(sum(A[b], j) - sum(A[b], i - 1))*p*q;
	}
	return ret;
/*
	ll ret = 0;
	for(int p=1;p<=n;p++)
		ret += 1LL*mu[p*b]*(n / p)*(m / p);
	return ret;
*/
}

bool tag[MAXN + 5];
void get(int k) {
	for(unsigned i=0;i<dv[k].size();i++)
		if( mu[dv[k][i]] ) B[++cnt] = dv[k][i], A[dv[k][i]] = cnt;
	for(int i=1;i<=cnt;i++) {
		tag[1] = true, smu[i][1] = 1;
		for(int j=2;j<=MAXN;j++) {
			if( !nprm[j] ) tag[j] = (B[i] % j);
			smu[i][j] = smu[i][j - 1] + (tag[j] ? mu[j] : 0);
			for(int p=1;j*prm[p]<=MAXN;p++) {
				nprm[j*prm[p]] = true, tag[j*prm[p]] = (tag[j] && tag[prm[p]]);
				if( j % prm[p] == 0 ) break;
			}
		}
		for(int j=1;j<=MAXN;j++) smu[i][j] *= mu[B[i]];
	}
}
int main() {
	int n, m, k;
	scanf("%d%d%d", &n, &m, &k);
	init(), get(k);
	
	ll ans = 0;
	for(unsigned i1=0;i1<dv[k].size();i1++) {
		int d = dv[k][i1]; if( !mu[d] ) continue;
		for(unsigned i2=0;i2<dv[d].size();i2++) {
			int b = dv[d][i2];
			for(unsigned i3=0;i3<dv[b].size();i3++) {
				int a = dv[b][i3];
				
				ans += mu[d] * mu[b / a] * solve(b, n / b, 1LL * m * a / b / d);
			}
		}
	}
	printf("%lld
", ans);
}

details

虽然 (summu) 不会溢出，但中间过程可能会溢出，所以还是要开 long long。