description
牛牛是一个热爱算法设计的高中生。在他设计的算法中，常常会使用带小数的数进行计算。牛牛认为，如果在 (k) 进制下，一个数的小数部分是纯循环的，那么它就是美的。
现在，牛牛想知道：对于已知的十进制数 (n) 和 (m)，在 (k) 进制下，有多少个数值上互不相等的纯循环小数，可以用分数 (frac{x}{y}) 表示，其中 (1leq x leq n, 1 leq y leq m)，且 (x, y) 是整数。
一个数是纯循环的，当且仅当其可以写成以下形式：
其中，(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) 也互质。因此答案为：
首先套路反演：
二度反演：
枚举四元组 ((a, b, d, k))，记 (n' = lfloorfrac{n}{b} floor, m'=lfloorfrac{m imes a}{d imes b} floor)，则我们相当于求：
后面部分很明显可以分块。考虑怎么快速求 (S_{b}(n) = sum_{i=1}^{n}mu(b imes i))，三度反演：
当 (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。