題目鏈接：循環之美
這道題感覺非常優美……能有一個這么優美的題面和較高的思維難度真的不容易……
為了表示方便，讓我先講一下兩個符號。\([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 循環之美