UOJ #221 【NOI2016】 循环之美

题目连接：循环之美

　　这道题感受很是优美……能有一个这么优美的题面和较高的思惟难度真的不容易……

　　为了表示方便，让我先讲一下两个符号。\([a]\)表示若是\(a\)为真，那么返回\(1\)，不然返回\(0\)； \(a \perp b\)表示\(a\)与\(b\)互质。

　　首先，咱们须要考虑一个分数要成为纯循环小数须要知足什么条件。

　　咱们先来回想一下，咱们是怎样使用除法来判断一个分数$\frac{x}{y}$是不是纯循环小数的。显然咱们是一路除下去，何时出现了相同的余数，那么这个数就是一个循环小数。若是第一个重复的余数是$x$，那么这个数就是纯循环小数。这种方法实际上就是存在一个数$l\neq 0$，使得：$$xk^l \equiv x (\bmod y)$$

　　又因为题目要求值不重复，那么咱们能够获得$x \perp y$。因此，咱们能够推出$k^l \equiv 1(\bmod y)$。因此咱们只需$k \perp y$便可。

　　因而，咱们其实是要求这个式子：

\begin{aligned}&\sum_{x=1}^{n} \sum_{y=1}^{m}[x\perp y][y \perp k]\\=&\sum_{y=1}^{m}[y\perp k]\sum_{x=1}^{n}[x\perp y]\end{aligned}

　　用一下莫比乌斯反演，能够获得：

\begin{aligned}&\sum_{y=1}^{m}[y\perp k] \sum_{x=1}^{n}\sum_{d|x,d|y}\mu(d)\\=&\sum_{y=1}^{m}[y\perp k]\sum_{d|y}^{n}\mu(d)\lfloor \frac{n}{d}\rfloor \end{aligned}

　　再经过换一下枚举顺序，可得：

\begin{aligned}&\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d} \rfloor \sum_{y=1}^{m}[d\ |\ y][y\perp k]\\=&\sum_{d=1}^{n}\mu(d)\lfloor \frac{n}{d}\rfloor \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor}[id\perp k]\\=& \sum_{d=1}^{n}[d\perp k]\mu(d)\lfloor \frac{n}{d}\rfloor \sum_{i=1}^{\lfloor \frac{m}{d} \rfloor}[i\perp k]\end{aligned}

　　最后一步用了互质的一个性质，那就是$[ab\perp c]=[a\perp c][b\perp c]$。这点应该很好理解吧。

　　咱们不妨设$f(n,k)=\sum_{i=1}^n[i\perp k]$，那么咱们只要在$O(1)$的时间内求出$f$，复杂度就是$O(n)$，就能够获得$84$分。实际上，咱们能够获得一个式子：

$$f(n,k)=\lfloor \frac{n}{k} \rfloor f(k,k)+f(n\bmod k,k)$$

　　这个应该也不难理解，由于$[a\perp b]=[a\bmod b \perp b]$

　　因而咱们就能够对$f(n,k)$预处理出$n\le k$的值，每次求值就变成$O(1)$了。因而$84$分到手了。

　　咱们考虑接下来该如何优化。因为$\lfloor \frac{m}{x} \rfloor$只有$\sqrt{m}$种取值，$\lfloor \frac{n}{x} \rfloor$只有$\sqrt{n}$种取值，因而咱们显然能够分段求和。而后，咱们就须要快速求出$\sum_{i=1}^n[i\perp k]\mu(i)$的值。

　　不妨设$g(n,k)=\sum_{i=1}^n[i\perp k]\mu(i)$，咱们来考虑一下这个函数如何快速求。咱们先考虑$k$的一个质因数$p$，那么$k$显然能够写成$p^cq$的形式。因为在$[1,n]$的范围中只有与$k$互质的才是有效值，那么若$x\perp k$，咱们能够获得$x\perp p$而且$x\perp q$。因而，咱们能够考虑从$x\perp q$的取值中减去$x$不与$p$互质的部分，就能够获得$x\perp k$的部分。这里若是不懂的话能够本身画一个$x\perp q$，$x\perp p$，$x\perp k$的关系图理解一下。

　　因为全部与$q$互质的数必定能够写成$p^xy(y\perp q)$的形式。那么咱们须要减去的数必定知足$x>0$。又因为当$x>1$时$\mu(p^xy)=0$，因此咱们只须要考虑$x=1$的状况便可。在这种状况下，咱们须要考虑的数就是$py(y\perp q)$的形式，因此咱们能够获得以下式子：\begin{aligned}  g(n,k)&=\sum_{i=1}^n[i\perp q]\mu(i)-\sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[py\perp q]\mu(py) \\&=g(n,q)- \sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp q]\mu(py)\end{aligned}

　　上面的最后一步是因为$p\perp q$，因此$py\perp q$只需在保证$y\perp q$便可。

　　咱们接着来考虑后一个式子。显然当$p\perp y$的时候$\mu(py)=\mu(p)\mu(y)$，不然$\mu(py)=0$。因而，咱们又获得了：

\begin{aligned} g(n,k)&=g(n,q)- \sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp p][y\perp q]\mu(p)\mu(y)\\&=g(n,q)-\mu(p)\sum_{y=1}^{\lfloor\frac{n}{p}\rfloor}[y\perp k]\mu(y)\\&=g(n,q)+g(\lfloor\frac{n}{p}\rfloor,k) \end{aligned}

　　因而咱们就能够递归求解了。容易发现边界状况就是$n=0$或者$k=1$。$n=0$的时候直接返回$0$就能够了，$k=1$的时候就是莫比乌斯函数的前缀和，用杜教筛求出来就能够了。因为每次递归要么$n$会变成$\lfloor \frac{n}{p} \rfloor$，有$O(\sqrt{n})$种取值；要么$p$少了一个质因数，有$\omega(k)$种取值，因此总共有$O(\omega(k)\sqrt{n})$种取值，记忆化搜索便可。其中$\omega(k)$表示$k$的不一样的质因子数目。因而最后的总复杂度为$O(\omega(k)\sqrt{n}+n^{\frac{2}{3}})$，能够经过此题。

　　两个细节：

　　能够在递归求$g$的时候不直接记$k$，而是记录$k$还剩下多少个质因数，方便存储；

　　记忆化的时候能够开哈希链表存储，用$map$也能够，或者干脆分$n,m$以及小于、大于根号的状况存储也能够。

　　第一次写这种题的题解，好像写得过于详细了(实际上是由于我什么都不会)，请大神们不要喷……

　　下面贴代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
#define maxn 5000010
#define mod 100007
#define maxk 2010
 
using namespace std;
typedef long long llg;
 
int n,m,k,pr_k[maxk],lp;
int s[maxn],ls,f[maxk];
int head[mod],next[maxn],to[maxn],tt;
int h1[11][mod],n1[maxn],t1[maxn],t2;
bool vis[maxn],huzhi[maxk];
llg ans,zto[maxn],zt1[maxn],mu[maxn];
 
int gcd(int a,int b){
    int r=a%b;
    while(r) a=b,b=r,r=a%b;
    return b;
}
 
void pre(){
    mu[1]=1;
    for(int i=2;i<maxn;i++){
        if(!vis[i]) s[++ls]=i,mu[i]=-1;
        for(int j=1;j<=ls && s[j]*i<maxn;j++){
            vis[s[j]*i]=1;
            if(i%s[j]) mu[i*s[j]]=-mu[i];
            else{mu[i*s[j]]=0;break;}
        }
    }
    for(int i=2;i<maxn;i++) mu[i]+=mu[i-1];
    for(int i=1;i<=k;i++){
        huzhi[i]=(gcd(i,k)==1);
        f[i]=f[i-1]+huzhi[i];
    }
}
 
llg suan(int x){return (x/k)*f[k]+f[x%k];}
llg solveu(int x){
    if(x<maxn) return mu[x];
    for(int i=head[x%mod];i;i=next[i]) if(to[i]==x) return zto[i];
    int now=++tt; to[tt]=x;next[tt]=head[x%mod];head[x%mod]=tt; zto[now]=1;
    for(int i=2,nt=0;nt<x;i=nt+1) nt=x/(x/i),zto[now]-=(nt-i+1)*solveu(x/i);
    return zto[now];
}
 
llg g(int x,int y){
    if(!x) return solveu(y);
    if(y<=1) return y;
    for(int i=h1[x][y%mod];i;i=n1[i]) if(t1[i]==y) return zt1[i];
    int now=++t2; t1[t2]=y;n1[t2]=h1[x][y%mod];h1[x][y%mod]=t2;
    return zt1[now]=g(x-1,y)+g(x,y/pr_k[x]);
}
 
int main(){
    File("a");
    scanf("%d %d %d",&n,&m,&k);
    pre();
    for(int i=1;s[i]<=k;i++) if(k%s[i]==0) pr_k[++lp]=s[i];
    for(int d=1,l=min(n,m),nt;d<=l;d=nt+1){
        nt=min(n/(n/d),m/(m/d));
        ans+=(g(lp,nt)-g(lp,d-1))*(llg)(n/d)*suan(m/d);
    }
    printf("%lld",ans);
    return 0;
}

　　 最后附一个BZOJ的提交网址：BZOJ 4652 循环之美