【UOJ #221】【NOI 2016】循环之美

http://uoj.ac/problem/221
因为(a)和(b)不互质时，(frac ab=frac{frac a{(a,b)}}{frac b{(a,b)}})，所以只用求(a)和(b)互质时的满足条件的个数。
(frac ab)在(k)进制下是纯循环小数，我们先假设循环节长度为(l)，这样(frac ab imes k^l-frac ab)的小数部分就是0，也就是这是个整数。
(frac{aleft(k^l-1 ight)}b)是个整数，就是说(b|aleft(k^l-1 ight))。又因为(aot b)，所以(b|left(k^l-1 ight))。
判断(frac ab)是否是(k)进制下的纯循环小数就转变成了判断是否存在一个(l,lgeq 0)，满足(left(k^l-1 ight)mod b=0)。
当(k)和(b)不互质时，对于所有的(l)，(b)和(k^l)有共同的质因子，又因为(k^lotleft(k^l-1 ight))，所以(b)和(k^l)的共同的质因子是(k^l-1)没有的，所以不存在(l)满足(b|left(k^l-1 ight))。
当(k)和(b)互质时，由欧拉定理：(k^{varphi(b)}mod b=1)，存在(l=varphi(b))。
这样答案就变成了：$$sum_{a=1}^nsum_{b=1}m[bot k][aot b]$$
化一波式子：

[egin{aligned} &sum_{a=1}^nsum_{b=1}^m[bot k][aot b]\ =&sum_{d=1}^nmu(d)sum_{a=1}^{leftlfloorfrac nd ight floor}sum_{b=1}^{leftlfloorfrac md ight floor}[bdot k]\ =&sum_{d=1}^nmu(d)[dot k]leftlfloorfrac nd ight floorsum_{b=1}^{leftlfloorfrac md ight floor}[bot k] end{aligned} ]

(sumlimits_{i=1}^n[iot k])很好预处理后(O(1))计算。
重点是怎么算(sumlimits_{d=1}^nmu(d)[dot k])？
先放宽限制，把([dot k])的限制去掉，这样就是对(Oleft(sqrt n+sqrt m ight))个下取整取值求(mu)的前缀和，可以先(Oleft(n^{frac 23} ight))大力杜教筛一波。
然后考虑(S(i,n))表示(1sim n)中与(k)的前(i)个质因子互质的数的(mu)值和，这样(S(i,n)=S(i-1,n)-muleft(p_i ight)Sleft(i,leftlfloorfrac n{p_i} ight floor ight))递推求出，时间复杂度(Oleft(left(sqrt n+sqrt m ight)log k ight))（假设(k)的不同质因子有(log k)个）。
最后对(Oleft(sqrt n+sqrt m ight))个不同的(leftlfloorfrac nd ight floor)和(leftlfloorfrac md ight floor)的取值进行分段求和即可。
时间复杂度(Oleft(n^{frac 23}+left(sqrt n+sqrt m ight)log k ight))。

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;

const int N = 1000003;

bool notp[N];
int mu[N], sum_mu[N], prime[N], Num = 0;

void Euler_shai() {
	mu[1] = sum_mu[1] = 1;
	for (int i = 2; i <= 1000000; ++i) {
		if (!notp[i]) prime[++Num] = i, mu[i] = -1;
		for (int j = 1; j <= Num && prime[j] * i <= 1000000; ++j) {
			notp[prime[j] * i] = true;
			if (i % prime[j] == 0) break;
			mu[prime[j] * i] = -mu[i];
		}
		sum_mu[i] = sum_mu[i - 1] + mu[i];
	}
}

int R[N], cnt = 0, sum_mu2[N], n, m, k, S[N], pr[N], prnum = 0;

const int mo = 2333333;

struct HashTable {
	int pos[mo], num[mo];
	
	void ins(int nu, int po) {
		int tmp = nu % mo;
		while (num[tmp]) {++tmp; if (tmp == mo) tmp = 0;}
		num[tmp] = nu;
		pos[tmp] = po;
	}
	
	int query(int nu) {
		int tmp = nu % mo;
		while (num[tmp] != nu) {++tmp; if (tmp == mo) tmp = 0;}
		return pos[tmp];
	}
} HT;

int sum[N], id[N];

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

void pre_sum() {
	for (int i = 1; i <= k; ++i) {
		sum[i] = sum[i - 1];
		if (gcd(i, k) == 1)
			++sum[i];
	}
}

int cal(int num) {return sum[k] * (num / k) + sum[num % k];}

int Sum(int num) {return num <= 1000000 ? S[id[num]] : S[HT.query(num)];}

int main() {
	scanf("%d%d%d", &n, &m, &k);
	Euler_shai();
	
	for (int i = 1, j = 1; i <= n && i <= m; ++i) {
		j = m / (m / i);
		i = n / (n / i);
		if (i > j) i = j;
		R[++cnt] = i;
		if (i > 1000000) HT.ins(i, cnt);
	}
	
	for (int i = 1; i <= cnt; ++i) {
		int num = R[i], &ret = S[i];
		if (num <= 1000000) {ret = sum_mu[num]; id[num] = i; continue;}
		ret = 1;
		for (int j = 2, pre = 1; j <= num; pre = j, ++j) {
			j = num / (num / j);
			ret -= 1ll * Sum(num / j) * (j - pre);
		}
	}
	
	for (int i = 1; i <= Num && prime[i] <= k; ++i)
		if (k % prime[i] == 0)
			pr[++prnum] = prime[i];
	
	for (int i = 1, pi = pr[1]; i <= prnum; pi = pr[++i])
		for (int j = 1; j <= cnt; ++j)
			S[j] += Sum(R[j] / pi);
	
	pre_sum();
	
	ll ans = 0;
	for (int tmp = 1, i = 1, j = 1; i <= n && i <= m; ++tmp, ++i) {
		j = m / (m / i);
		i = n / (n / i);
		if (i > j) i = j;
		ans += 1ll * (S[tmp] - S[tmp - 1]) * (n / i) * cal(m / i);
	}
	
	printf("%lld
", ans);
	return 0;
}