【NOI2016】循环之美，mobius反演+杜教筛

思路：
对于xy,我们肯定让gcd(x,y)=1，因为这样计算不会重复，而且比不互质时更优，且只要余数出现重复就是有循环节了，根据十进制小数转k进制的方法，我们可以得到
x⋅ks=x(mody)⇒ks=1(mody)或y|x
因为gcd(x,y)=1，所以gcd(ks,y)=1，也就是gcd(k,y)=1
我们要求的式子就是∑i=1n∑j=1m[gcd(i,j)=1][gcd(j,k)=1]
把[gcd(i,j)=1]移到后面，反演一下，再在前面更换指标
∑dμ(d)⌊nd⌋∑j=1⌊md⌋[gcd(dj,k)=1]
(然后我就反演了半天，最后被reflash和mrazer一眼秒了)
因为gcd(dj,k)=1⇔gcd(d,k)=1,gcd(j,k)=1，所以
∑gcd(d,k)=1μ(d)⌊nd⌋F(⌊md⌋)，其中F(m)=∑i=1m[gcd(i,k)=1]
关于F函数，因为gcd(i,k)=gcd(i+k,k)，所以很多都是重复的，可以预处理k以内的F，然后就能O(1)算了
（我原本没预处理，暴力算的，在uoj跑76分，问了reflash才知道可以O(1)算的）
前面部分根据d的取值分块算，现在问题就在于如何快速计算∑gcd(d,k)=1μ(d)
%一发reflash
(这个我真没想到，卡在这里挺久的，因为没有把μ(di)拆开看）
设g(n)=∑i=1n[gcd(i,k)=1]μ(i),M(n)=∑i=1nμ(i)
则g(n)=M(n)−∑d=2n∑i=1n[gcd(i,k)=d]μ(i)=M(n)−∑d=2,d|kn∑i=1⌊nd⌋[gcd(i,kd)=1]μ(di)=M(n)−∑d=2,d|kn∑i=1⌊nd⌋[gcd(i,kd)=1]μ(d)μ(i)[gcd(i,d)=1]
两个互质条件合并一下，就是[gcd(i,k)=1]
g(n)=M(n)−∑d=2,d|knμ(d)g(⌊nd⌋)
非常像杜教筛的处理方法，记忆化即可
小范围的时候可以暴力
代码：

#include<cstdio>#include<iostream>#define LL long long#define ui unsigned intusing namespace std;int gcd(int x,int y){return y?gcd(y,x%y):x;}int n,m,k;const int lim=2000000;int prime[lim/10+5],f[2005],py[2005];ui pr[2005];short mu[lim+5],sum[lim+5];char vis[lim+5];void init(){    mu[1]=1;    int t,i,j;    for (i=2;i<=lim;++i)    {        if (!vis[i])            prime[++prime[0]]=i,            mu[i]=-1;        for (j=1;j<=prime[0]&&(t=i*prime[j])<=lim;++j)        {            vis[t]=1;            if (i%prime[j]) mu[t]=-mu[i];            else break;         }    }    for (int i=1;i<=lim;++i) sum[i]=sum[i-1]+mu[i];    for (int i=2;i<=k;++i)        if (k%i==0&&mu[i]) pr[++pr[0]]=i;    for (int i=1;i<=k;++i)    {        f[i]=f[i-1];py[i]=py[i-1];        if (gcd(i,k)==1) ++f[i],py[i]+=mu[i];    }}LL S(ui m){return 1LL*(m/k)*f[k]+f[m%k];}const int P=999983;struct node{    int tot,first[P+5],next[200005],ha[200005];    ui data[200005];    void push(ui n,int t)    {        int p=n%P;        ha[++tot]=t;        data[tot]=n;        next[tot]=first[p];        first[p]=tot;    }    int get(ui n)    {        for (int i=first[n%P];i;i=next[i])            if (data[i]==n) return i;        return 0;    }}mp_mu,mp;int cal_mu(ui n){    if (n<=lim) return sum[n];    int id=mp_mu.get(n);    if (id) return mp_mu.ha[id];    int t=1;    for (ui last,p,i=2;i<=n;i=last+1)        last=n/(p=n/i),        t-=(last-i+1)*cal_mu(p);    mp_mu.push(n,t);    return t;}int cal(ui n){    if (n<=k) return py[n];    int id=mp.get(n);    if (id) return mp.ha[id];    int t=cal_mu(n);    for (int i=1;i<=pr[0]&&pr[i]<=n;++i)        t-=mu[pr[i]]*cal(n/pr[i]);    mp.push(n,t);    return t;}main(){    scanf("%d%d%d",&n,&m,&k);    init();    LL ans=0;    int ls=0,rs,ll=min(n,m);    for (ui p,q,last,i=1;i<=ll;i=last+1)        last=min((ui)n/(p=n/i),(ui)m/(q=m/i)),        ans+=((rs=cal(last))-ls)*S(q)*p,        ls=rs;    printf("%lld\n",ans);}