LibreOJ 2085 「NOI2016」循环之美

首先，为了去重，我们可以增加一个限制：必须是最简分数。
然后我们看看什么情况会产生纯循环小数。

在 10 进制下，当 gcd(x,y)=1,gcd(x,10)=1 时，yx 是纯循环小数。
于是我们猜测当 gcd(x,y)=1,gcd(x,k)=1 时，yx 在 k 进制下是纯循环小数。

证明：
yx 在 k 进制下是纯循环小数，当且仅当 ∃l≠0,yx−⌊yx⌋=yklx−⌊yklx⌋。
这玩意看着很像取模的形式，所以两边同乘 x，变成 y≡ykl(modx)。

由于前面说了 gcd(x,y)=1，所以 kl≡1(modx)。
于是由欧拉定理得 gcd(k,x)=1，一个特解为 l=φ(x)。

证毕。

问题变成了 n∑i=1m∑j=1[gcd(i,j)=1][gcd(j,k)=1]。

n∑i=1m∑j=1[gcd(i,j)=1][gcd(j,k)=1]=m∑j=1[gcd(j,k)=1]n∑i=1[gcd(i,j)=1]=m∑j=1[gcd(j,k)=1]n∑i=1∑d|i,d|jμ(d)=m∑j=1[gcd(j,k)=1]∑d|j⌊nd⌋μ(d)=min(n,m)∑d=1⌊nd⌋μ(d)∑jd≤m[gcd(jd,k)=1]=min(n,m)∑d=1⌊nd⌋[gcd(d,k)=1]μ(d)∑jd≤m[gcd(j,k)=1]

设 f(n)=n∑i=1[gcd(i,k)=1]。
由于 gcd(i,k)=gcd(i+k,k)，有 f(n)=⌊nk⌋φ(k)+f(nmodk)。
发现 k 非常小，暴力求出 x∈[0,k] 的 f(x) 值即可。

于是有 min(n,m)∑d=1⌊nd⌋[gcd(d,k)=1]μ(d)∑jd≤m[gcd(j,k)=1]=min(n,m)∑d=1⌊nd⌋f(⌊md⌋)[gcd(d,k)=1]μ(d)

如果数论分块的话，就剩下后面的两项。
考虑怎么算它们的前缀和。

设 S(n,k)=n∑i=1[gcd(i,k)=1]μ(i)。
有 S(n,k)=n∑i=1[gcd(i,k)=1]μ(i)=n∑i=1μ(i)∑d|i,d|kμ(d)=∑d|kμ(d)∑id≤nμ(id)

乍一看推不下去了。
但是我们反观 μ 函数，根据定义容易发现当 gcd(x,y)≠1 时，μ(xy)=0。

所以有 S(n,k)=∑d|kμ(d)∑id≤nμ(id)=∑d|kμ(d)∑id≤n[gcd(i,d)=1]μ(id)=∑d|kμ2(d)∑id≤n[gcd(i,d)=1]μ(i)=∑d|kμ2(d)S(⌊nd⌋,d)

然后考虑递归求解 S，容易发现 S(n,1)=n∑i=1μ(i)，杜教筛即可。

于是最后对于 min(n,m)∑d=1⌊nd⌋f(⌊md⌋)[gcd(d,k)=1]μ(d) 数论分块即可。

略微卡常，少开 long long，算 S 的时候如果系数是 0 就不用往下算。
unordered_map 居然不能用 pair，偷懒换成了 map（

代码：

#include <cstdio>
#include <utility>
#include <algorithm>
#include <map>
#include <tr1/unordered_map>
using namespace std;
using namespace tr1;
const int N = 1e9;
const int K = 2e3;
const int CNT = 1e7;
int n,m,k;
int vis[CNT + 5],prime[CNT + 5],cnt,mu[CNT + 5];
int f[K + 5];
unordered_map<int,int> w;
map<pair<int,int>,long long> mem;
long long ans;
int calc(int n)
{
    if(n <= CNT)
        return mu[n];
    if(w.count(n))
        return w[n];
    long long ret = 1;
    for(register int l = 2,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret -= (long long)(r - l + 1) * calc(n / l);
    }
    return w[n] = ret;
}
int F(int n)
{
    return f[k] * (n / k) + f[n % k];
}
long long S(int n,int k)
{
    if(n == 0 || n == 1)
        return n;
    if(k == 1)
        return calc(n);
    if(mem.count(make_pair(n,k)))
        return mem[make_pair(n,k)];
    long long ret = 0;
    for(int i = 1;i * i <= k;++i)
        if(!(k % i))
        {
            if(mu[i] - mu[i - 1])
                ret += (mu[i] - mu[i - 1]) * (mu[i] - mu[i - 1]) * S(n / i,i);
            if(i ^ (k / i) && mu[k / i] - mu[k / i - 1])
                ret += (mu[k / i] - mu[k / i - 1]) * (mu[k / i] - mu[k / i - 1]) * S(n / (k / i),k / i);
        }
    return mem[make_pair(n,k)] = ret;
}
int main()
{
    mu[1] = 1;
    for(register int i = 2;i <= CNT;++i)
    {
        if(!vis[i])
            mu[prime[++cnt] = i] = -1;
        for(register int j = 1;j <= cnt && i * prime[j] <= CNT;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
                break;
            else
                mu[i * prime[j]] = -mu[i];
        }
        mu[i] += mu[i - 1];
    }
    scanf("%d%d%d",&n,&m,&k);
    for(register int i = 1;i <= k;++i)
        f[i] = f[i - 1] + (__gcd(i,k) == 1);
    for(register int l = 1,r;l <= min(n,m);l = r + 1)
    {
        r = min(n / (n / l),m / (m / l));
        ans += (S(r,k) - S(l - 1,k)) * (n / l) * F(m / l);
    }
    printf("%lld\n",ans);
}