BZOJ 4652: [Noi2016]循环之美


BZOJ提交链接:
http://www.lydsy.com/JudgeOnline/problem.php?id=4652

结论:一个既约分数a/b在k进制下是纯循环小数的充要条件是gcd(b,k)=1

先推一波式子
$$\sum_{i=1}^{n} \sum_{j=1}^{m} [gcd(i,j)=1][gcd(j,k)=1]$$
$$=\sum_{j=1}^{m} [gcd(j,k)=1] \sum_{i=1}^{n} [gcd(i,j)=1]$$
$$=\sum_{j=1}^{m} [gcd(j,k)=1] \sum_{i=1}^{n} \sum_{d|i,d|j} \mu(d)$$
$$=\sum_{j=1}^{m} [gcd(j,k)=1] \sum_{d|j} \mu(d) \lfloor \frac{n}{d} \rfloor$$
$$=\sum_{j=1}^{m} \sum_{d|j} \lfloor \frac{n}{d} \rfloor [gcd(j,k)=1] \mu(d)$$
$$=\sum_{d=1}^{m} \mu(d) \lfloor \frac{n}{d} \rfloor \sum_{d|j,j<=m} [gcd(j,k)=1]$$
$$=\sum_{d=1}^{m} \mu(d) \lfloor \frac{n}{d} \rfloor \sum_{x=1}^{\lfloor \frac{m}{d} \rfloor} [gcd(xd,k)=1]$$
$$=\sum_{d=1}^{m} [gcd(d,k)=1] \mu(d) \lfloor \frac{n}{d} \rfloor \sum_{x=1}^{\lfloor \frac{m}{d} \rfloor} [gcd(x,k)=1]$$

以上一些步骤用到了$$[gcd(xy,k)=1]=[gcd(x,k)=1][gcd(y,k)=1]$$的性质

如何计算这个式子呢…

考虑如何计算 $$\sum_{k=1}^a [gcd(k, b)=1]$$

我们用 $$s(a,b)$$ 表示 $$\sum_{k=1}^a [gcd(k, b)=1]$$ 的值
我们把 k 按除以b的余数分类。
那么
$$s(a,b)=\lfloor\frac{a}{b}\rfloor s(b, b)+s(a\%b,b)$$
k的值很小,我们可以预处理一些a小于b的s(a,b)。
这样就可以做到O(n)计算了

继续优化

回忆一下推出的式子
$$\sum_{d=1}^{m} [gcd(d,k)=1] \mu(d) \lfloor \frac{n}{d} \rfloor \sum_{x=1}^{\lfloor \frac{m}{d} \rfloor} [gcd(x,k)=1]$$
可以写成
$$\sum_{d=1}^{m} [gcd(d,k)=1] \mu(d) \lfloor \frac{n}{d} \rfloor s(\lfloor \frac{m}{d} \rfloor,k)$$
$$\lfloor \frac{n}{d} \rfloor$$ 和 $$\lfloor \frac{m}{d} \rfloor$$ 都取值很少,可以数论分块。
那么问题变成,如何快速计算$$\sum_{d=1}^{m} [gcd(d,k)=1] \mu(d)$$
考虑杜教筛。设
$$f(n)=[gcd(n,k)=1] \mu(n)$$
$$F[n]=\sum_{i=1}^n f(i)$$
$$g(n)=[gcd(n,k)=1]$$
可以得到
$$(f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})=\sum_{d|n}[gcd(d,k)=1]\mu(d)[gcd(\frac{n}{d},k)=1]=[gcd(n,k)=1]\sum_{d|n}\mu(d)=[gcd(n,k)=1][n=1]=[n=1]$$

$$\sum_{i=1}^n (f*g)(i)=1$$

$$\sum_{i=1}^n \sum_{d|i} f(\frac{i}{d})g(d)=1$$
$$\sum_{d=1}^n g(d) \sum_{d|i,i<=n} f(\frac{i}{d})=1$$
$$\sum_{d=1}^n [gcd(d,k)=1] F(\lfloor \frac{n}{d} \rfloor)=1$$
$$F(n)=1-\sum_{d=2}^n [gcd(d,k)=1] F(\lfloor \frac{n}{d} \rfloor)$$
再进行数论分块,结合前面提到的计算s(a,b)的方法,就能够计算出答案。
下面附上AC代码

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cstdio>
#include <map>

using namespace std;

typedef long long ll;

int n, m, k;

const int maxk = 2010;
const int MAX = 1000000;

int s_[maxk][maxk];
int f[MAX+10], F[MAX+10];
int isnp[MAX+10], p[MAX+10];
int tot = 0;

map<int, int> mp;

ll s(int a, int b) {
    return s_[a%b][b]+(a/b)*s_[b][b];
}
int gcd(int a, int b) {
    return !b?a:gcd(b, a%b);
}
void init() {
    isnp[1] = 1; f[1] = 1;
    for (int i = 2; i <= MAX; i++) {
        if (!isnp[i]) {
            p[++tot] = i;
            f[i] = -(k%i!=0);
        }
        for (int j = 1; j <= tot && p[j]*i <= MAX; j++) {
            isnp[i*p[j]] = 1;
            if (i % p[j]) {
                f[i*p[j]] = f[i]*f[p[j]];
            } else {
                f[i*p[j]] = 0;
                break;
            }
        }
    }
    for (int i = 1; i <= MAX; i++) F[i] = F[i-1]+f[i];
}
int getF(int x) {
    if (x <= MAX) return F[x];
    if (mp.find(x) != mp.end()) return mp[x];
    int ret = 1;
    for (int d = 2; d <= x;) {
        int nxt = x/(x/d);
        ret -= getF(x/d)*(s(nxt, k)-s(d-1,k));
        d = nxt+1;
    }
    return mp[x]=ret;
}

int main() {
    scanf("%d%d%d", &n, &m, &k);
    for (int j = 1; j <= k; j++) {
        for (int i = 1; i <= k; i++) {
            s_[i][j] = s_[i-1][j]+(gcd(i, j)==1);
        }
    }
    init();
    ll ret = 0;
    for (int d = 1; d <= m;) {
        int nxt1 = m, nxt2 = m;
        if (n/d) nxt1 = n/(n/d);
        if (m/d) nxt2 = m/(m/d);
        int nxt = min(nxt1, nxt2);
        ret += 1ll*(getF(nxt)-getF(d-1))*(n/d)*s(m/d, k);
        d = nxt+1;
    }
    cout << ret << endl;
    return 0;
}