[NOI2016]循环之美(杜教筛)

首先要求每个数互不相等，故有$x\perp y$。

可以发现$\frac{x}{y}$在$k$进制下为纯循环小数的充要条件为$x\cdot k^{len}\equiv x(mod\ y)$，即$y\perp k$。

接下来进行经典的推导：
$$\begin{aligned}&\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[i\perp j][j\perp k]\\=&\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum_{d|i,d|j}\mu(d)[j\perp k]\\=&\sum\limits_{d=1}^{min(n,m)}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}[jd\perp k]\\=&\sum\limits_{d=1}^{min(n,m)}\mu(d)\lfloor\frac{n}{d}\rfloor\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}[j\perp k][d\perp k]\end{aligned}$$

令$f(n)=\sum\limits_{i=1}^{n}[i\perp k]$，由于$gcd(i,k)=gcd(i-k,k)$，故$f(n)=\lfloor\frac{n}{k}\rfloor\cdot f(k)+f(n\%k)$

再令$g(i)=\mu(i)[i\perp k]$，则答案为$\sum\limits_{d=1}^{min(n,m)}g(d)f(\lfloor\frac{m}{d}\rfloor)\lfloor\frac{n}{d}\rfloor$

令$G(n,k)$为$g()$的前缀和，同样进行根号优化：
$$G(n,k)=\sum\limits_{i=1}^{n}\mu(i)[i\perp k]=\sum\limits_{i=1}^{n}\mu(i)\sum\limits_{d|i,d|k}\mu(d)=\sum\limits_{d|k}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(id)$$
注意到$\mu(id)=[i\perp d]\ ?\ 0:\mu(i)\cdot \mu(d)$，故
$$G(n,k)=\sum\limits_{d|k}\mu^2(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)[i\perp d]=\sum\limits_{d|k}\mu^2(d)G(\lfloor\frac{n}{d}\rfloor,d)$$

$G$函数可以通过递归记忆化求出。由于$G(a,b)$中$a$最多有$\sqrt{n}$种取值，$b$最多有$\sigma_0(k)$种取值，每次转移的复杂度也是$\sigma_0(k)$的，故总复杂度为$O(n^{\frac{2}{3}}+\sigma_0^2(k))$，事实上远远达不到这个上界。

在做除$\mu$和$\varphi$之外的杜教筛时，用map会简洁得多，有时候（可能询问到的n不能确定时）必须用map。

此题还有一种更优类似洲阁筛的做法，能处理$k\leqslant 10^{18}$的问题，复杂度为$O(n^{\frac{2}{3}}+\omega(k)\sqrt{n})$。

https://blog.sengxian.com/solutions/bzoj-4652

#include<map>
#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
typedef long long ll;
using namespace std;

const int N=5000010,K=2010;
ll ans;
bool b[N];
int tot,w,n,m,k,F[K],p[N],miu[N];
map<int,int>G[K],Miu;

int gcd(int a,int b){ return b ? gcd(b,a%b) : a; }
int f(int n){ return (n/k)*F[k]+F[n%k]; }

void pre(int n){
    miu[1]=1;
    rep(i,2,n){
        if (!b[i]) p[++tot]=i,miu[i]=-1;
        for (int j=1; j<=tot && i*p[j]<=n; j++){
            b[i*p[j]]=1;
            if (i%p[j]==0) { miu[i*p[j]]=0; break; }
                else miu[i*p[j]]=-miu[i];
        }
    }
    rep(i,2,n) miu[i]+=miu[i-1];
}

int Mu(int n){
    int res=1;
    if (n<=w) return miu[n];
    if (Miu.count(n)) return Miu[n];
    for (int i=2,lst; i<=n; i=lst+1)
        lst=n/(n/i),res-=Mu(n/i)*(lst-i+1);
    return Miu[n]=res;
}

int g(int n,int k){
    int res=0;
    if (!n) return 0;
    if (G[k].count(n)) return G[k][n];
    if (k==1) return Mu(n);
    for (int d=1; d*d<=k; d++)
        if (k%d==0){
            if (miu[d]-miu[d-1]) res+=g(n/d,d);
            if (d*d==k) continue;
            int t=k/d;
            if (miu[t]-miu[t-1]) res+=g(n/t,t);
        }
    return G[k][n]=res;
}

int main(){
    freopen("cycle.in","r",stdin);
    freopen("cycle.out","w",stdout);
    scanf("%d%d%d",&n,&m,&k); w=min(max(n,m),5000000); pre(w);
    rep(i,1,k) F[i]=F[i-1]+(gcd(i,k)==1);
    for (int i=1,lst; i<=m && i<=n; i=lst+1)
        lst=min(n/(n/i),m/(m/i)),ans+=1ll*(g(lst,k)-g(i-1,k))*f(m/i)*(n/i);
    printf("%lld\n",ans);
    return 0;
}