首先要求每个数互不相等，故有$xperp y$。
可以发现$frac{x}{y}$在$k$进制下为纯循环小数的充要条件为$xcdot k^{len}equiv x(mod y)$，即$yperp k$。
接下来进行经典的推导：
$$egin{aligned}&sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[iperp j][jperp k]\=&sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sum_{d|i,d|j}mu(d)[jperp k]\=&sumlimits_{d=1}^{min(n,m)}mu(d)sumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor}[jdperp k]\=&sumlimits_{d=1}^{min(n,m)}mu(d)lfloorfrac{n}{d}
floorsumlimits_{j=1}^{lfloorfrac{m}{d}
floor}[jperp k][dperp k]end{aligned}$$
令$f(n)=sumlimits_{i=1}^{n}[iperp k]$，由于$gcd(i,k)=gcd(i-k,k)$，故$f(n)=lfloorfrac{n}{k} floorcdot f(k)+f(n\%k)$
再令$g(i)=mu(i)[iperp k]$，则答案为$sumlimits_{d=1}^{min(n,m)}g(d)f(lfloorfrac{m}{d} floor)lfloorfrac{n}{d} floor$
令$G(n,k)$为$g()$的前缀和，同样进行根号优化：
$$G(n,k)=sumlimits_{i=1}^{n}mu(i)[iperp k]=sumlimits_{i=1}^{n}mu(i)sumlimits_{d|i,d|k}mu(d)=sumlimits_{d|k}mu(d)sumlimits_{i=1}^{lfloorfrac{n}{d}
floor}mu(id)$$
注意到$mu(id)=[iperp d] ? 0:mu(i)cdot mu(d)$，故
$$G(n,k)=sumlimits_{d|k}mu^2(d)sumlimits_{i=1}^{lfloorfrac{n}{d}
floor}mu(i)[iperp d]=sumlimits_{d|k}mu^2(d)G(lfloorfrac{n}{d}
floor,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。
此题还有一种更优类似洲阁筛的做法，能处理$kleqslant 10^{18}$的问题，复杂度为$O(n^{frac{2}{3}}+omega(k)sqrt{n})$。
https://blog.sengxian.com/solutions/bzoj-4652
1 #include<map> 2 #include<cstdio> 3 #include<algorithm> 4 #define rep(i,l,r) for (int i=(l); i<=(r); i++) 5 typedef long long ll; 6 using namespace std; 7 8 const int N=5000010,K=2010; 9 ll ans; 10 bool b[N]; 11 int tot,w,n,m,k,F[K],p[N],miu[N]; 12 map<int,int>G[K],Miu; 13 14 int gcd(int a,int b){ return b ? gcd(b,a%b) : a; } 15 int f(int n){ return (n/k)*F[k]+F[n%k]; } 16 17 void pre(int n){ 18 miu[1]=1; 19 rep(i,2,n){ 20 if (!b[i]) p[++tot]=i,miu[i]=-1; 21 for (int j=1; j<=tot && i*p[j]<=n; j++){ 22 b[i*p[j]]=1; 23 if (i%p[j]==0) { miu[i*p[j]]=0; break; } 24 else miu[i*p[j]]=-miu[i]; 25 } 26 } 27 rep(i,2,n) miu[i]+=miu[i-1]; 28 } 29 30 int Mu(int n){ 31 int res=1; 32 if (n<=w) return miu[n]; 33 if (Miu.count(n)) return Miu[n]; 34 for (int i=2,lst; i<=n; i=lst+1) 35 lst=n/(n/i),res-=Mu(n/i)*(lst-i+1); 36 return Miu[n]=res; 37 } 38 39 int g(int n,int k){ 40 int res=0; 41 if (!n) return 0; 42 if (G[k].count(n)) return G[k][n]; 43 if (k==1) return Mu(n); 44 for (int d=1; d*d<=k; d++) 45 if (k%d==0){ 46 if (miu[d]-miu[d-1]) res+=g(n/d,d); 47 if (d*d==k) continue; 48 int t=k/d; 49 if (miu[t]-miu[t-1]) res+=g(n/t,t); 50 } 51 return G[k][n]=res; 52 } 53 54 int main(){ 55 freopen("cycle.in","r",stdin); 56 freopen("cycle.out","w",stdout); 57 scanf("%d%d%d",&n,&m,&k); w=min(max(n,m),5000000); pre(w); 58 rep(i,1,k) F[i]=F[i-1]+(gcd(i,k)==1); 59 for (int i=1,lst; i<=m && i<=n; i=lst+1) 60 lst=min(n/(n/i),m/(m/i)),ans+=1ll*(g(lst,k)-g(i-1,k))*f(m/i)*(n/i); 61 printf("%lld ",ans); 62 return 0; 63 }