【UOJ221】【BZOJ4652】【NOI2016】循环之美（莫比乌斯反演，杜教筛）

Description

牛牛是一个热爱算法设计的高中生。在他设计的算法中，常常会使用带小数的数进行计算。牛牛认为，如果在 k 进制下，一个数的小数部分是纯循环的，那么它就是美的。现在，牛牛想知道：对于已知的十进制数 n 和 m，在 kk 进制下，有多少个数值上互不相等的纯循环小数，可以用分数 xy 表示,其中 1≤x≤n,1≤y≤m，且 x,y是整数。一个数是纯循环的，当且仅当其可以写成以下形式：a.c1˙c2c3…cp-1cp˙其中，a 是一个整数，p≥1；对于 1≤i≤p，ci是 kk 进制下的一位数字。例如，在十进制下，0.45454545……=0.4˙5˙是纯循环的，它可以用 5/11、10/22 等分数表示；在十进制下，0.1666666……=0.16˙则不是纯循环的，它可以用 1/6 等分数表示。需要特别注意的是，我们认为一个整数是纯循环的，因为它的小数部分可以表示成 0 的循环或是 k-1 的循环；而一个小数部分非 0 的有限小数不是纯循环的。

Solution

感觉这题的题面比较巧妙，至少一眼看不出是杜教筛。。。

有几个地方推不动了，果然太菜了，，不看题解根本做不出来啊。。。

XJB猜结论发现题目其实就是求：

$$\sum_{i=1}^{n}\sum_{j=1}^{m}[(i,j)=1][(j,k)=1]$$
然后颓柿子：
$$=\sum_{j=1}^m[(j,k)=1]\sum_{i=1}^n [(i,j)=1]$$
$$=\sum_{j=1}^m [(j,k)=1]\sum_{i=1}^n \sum_{d|(i,j)}\mu(d)$$
$$=\sum_{d=1}^n \mu(d)\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}[(jd,k)=1]\lfloor\frac{n}{d}\rfloor$$
$$=\sum_{d=1}^n [(d,k)=1]\mu(d)\lfloor\frac{n}{d}\rfloor\sum_{j=1}^{\lfloor \frac{m}{d}\rfloor}[(j,k)=1]$$
设$g(x)=\sum_{i=1}^x [(i,k)=1]$，那么有$g(x)=\lfloor\frac{x}{k}\rfloor g(k)+g(x\bmod k)$
有了这个后，继续颓：
$$=\sum_{d=1}^n [(d,k)=1]\mu(d)\lfloor\frac{n}{d}\rfloor g(\lfloor \frac{m}{d}\rfloor)$$

然后我就不会做了。。。以前没怎么遇到过这种递归解决的题目，真的挺巧(sang4)妙
(bing4)的。

我们设$f(n,k)=\sum_{d=1}^n [(d,k)=1]\mu(d)$，$k=p^kq$，$p$为质数且$(p,q)=1$，那么有：

$$f(n,k)=f(n,q)-\sum_{d=1}^{\lfloor\frac{n}{p}\rfloor}[(dp,q)=1]\mu(dp)$$
（减去与q互质单不和p互质的数的贡献）

$$f(n,k)=f(n,q)-\sum_{d=1}^{\lfloor\frac{n}{p}\rfloor}[(d,q)=1]\mu(dp)$$
（$(p,q)=1$）

$$=f(n,q)+\sum_{d=1}^{\lfloor\frac{n}{p}\rfloor}[(d,k)=1]\mu(d)$$

（计入贡献$d$一定和$p$互质（因为如果不互质则$\mu(dp)=0$）。$\mu(dp)$中除掉一个$p$后反号）

即：

$$f(n,k)=f(n,q)-f(\lfloor\frac{n}{p}\rfloor,k)$$

直接整除分块并递归求解，其中$\mu$前缀和用杜教筛即可。

Source

/**************************************
 * Au: Hany01
 * Prob: [BZOJ4652][NOI2016] 循环之美
 * Date: Mar 29th, 2018
 * Email: hany01@foxmail.com
**************************************/

#include<bits/stdc++.h>
#include<ext/pb_ds/hash_policy.hpp>
#include<ext/pb_ds/assoc_container.hpp>

using namespace std;
using namespace __gnu_pbds;

typedef long long LL;
typedef pair<int, int> PII;
typedef vector<int> VI;
#define File(a) freopen(a".in", "r", stdin), freopen(a".out", "w", stdout)
#define rep(i, j) for (register int i = 0, i##_end_ = j; i < i##_end_; ++ i)
#define For(i, j ,k) for (register int i = (j), i##_end_ = (k); i <= i##_end_; ++ i)
#define Fordown(i, j, k) for (register int i = (j), i##_end_ = (k); i >= i##_end_; -- i)
#define Set(a, b) memset(a, b, sizeof(a))
#define SZ(a) ((int)(a.size()))
#define ALL(a) a.begin(), a.end()
#define pb(a) push_back(a)
#define mp(a, b) make_pair(a, b)
#define INF (0x3f3f3f3f)
#define INF1 (2139062143)
#define Mod (1000000007)
#define y1 wozenmezhemecaia 
#ifdef hany01
#define debug(...) fprintf(stderr, __VA_ARGS__)
#else
#define debug(...)
#endif

template<typename T> inline bool chkmax(T &a, T b) { return a < b ? a = b, 1 : 0; }
template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }

inline int read() {
    register char c_; register int _, __;
    for (_ = 0, __ = 1, c_ = getchar(); !isdigit(c_); c_ = getchar()) if (c_ == '-')  __ = -1;
    for ( ; isdigit(c_); c_ = getchar()) _ = (_ << 1) + (_ << 3) + (c_ ^ 48);
    return _ * __;
}

const int maxn = 100000005, maxk = 2005, maxN = 100005;

int n, m, k, mu[maxN], pr[maxN], cnt, np[maxN], md[maxN], G[maxk], N, md1[maxN];
LL Ans;
cc_hash_table<LL, LL> ht;
cc_hash_table<int, LL> ht_;

inline void getmu()
{
    mu[1] = 1;
    For(i, 2, N) {
        if (!np[i])
            pr[++ cnt] = i, mu[i] = -1, md[i] = i, md1[i] = 1;
        for (register int j = 1; j <= cnt && pr[j] * i <= N; ++ j) {
            np[i * pr[j]] = 1, md[i * pr[j]] = pr[j];
            if (!(i % pr[j])) { mu[i * pr[j]] = 0, md1[i * pr[j]] = md1[i]; break; }
            mu[i * pr[j]] = -mu[i], md1[i * pr[j]] = i;
        }
    }
    For(i, 2, N) mu[i] += mu[i - 1];
}

inline void getg()
{
    For(i, 1, k) G[i] = G[i - 1] + (__gcd(i, k) == 1);
}

inline LL g(int x) { return x / k * G[k] + G[x % k]; }

LL S(int n)
{
    if (n < N) return mu[n];
    if (ht_[n]) return ht_[n];
    LL Ans = 1;
    for (int l = 2, r; l <= n; l = r + 1)
        r = n / (n / l), Ans -= (LL)(r - l + 1) * S(n / l);
    return ht_[n] = Ans;
}

LL f(int n, int k)
{
    if (!n) return 0;
    if (k == 1) return S(n);
    register LL tmp = (LL)k * (LL)maxn + (LL)n;
    if (!ht[tmp]) ht[tmp] = f(n, md1[k]) + f(n / md[k], k);
    return ht[tmp];
}

int main()
{
#ifdef hany01
    File("bzoj4652");
#endif

    //Input
    n = read(), m = read(), k = read();

    N = 100000;
    getmu();

    getg();

    int nm = min(n, m);
    for (int l = 1, r; l <= nm; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        Ans += (f(r, k) - f(l - 1, k)) * (LL)(n / l) * g(m / l);
    }

    cout << Ans << endl;

    return 0;
}
//坡谓西湖，正如西子，浓抹淡妆临镜台。
//    -- 刘过《六州歌头·寄稼轩承旨》