原文链接：[sol] P1587 [NOI2016] 循环之美 - icyM3tra's Blog
此为备份。
0. 前言
提供一种式子较为简洁且不算难写的做法。
以及介绍一种较为优秀的 dp 版杜教筛实现。
1. 式子推导
我们知道纯循环小数一定可以写成分母形如 \(k^a-1\) 的分数。
因此约分后的分母一定是 \(k^a-1\) 的约数。那么答案即为：
方括号中的两个条件分别是 \(x,y\) 互质、\(y|k^a-1\) 有正整数解 \(a\)。
容易证明后一条等价于 \(y\perp k\)。
考虑对其进行莫比乌斯反演：
注意到可以对 \(d\) 做整除分块。
接下来我们希望实现：
- 对 \(f(n)=\mu(n)[n\perp k]\) 求和；
- 对 \(g(n)=[n\perp k]\) 求和。
后者比较好处理。
容易发现 \(k\) 是 \(g\) 的周期。
我们 \(O(k\log k)\) 暴力计算 \(g(1)\sim g(k)\)，并预处理这一段前缀和，然后即可 \(O(1)\) 实现对 \(g\) 求和。
前者事实上是一个积性函数求和问题。
恰好我们注意到 \(f*g=\varepsilon\)，因此可以直接用杜教筛处理。
注意不到的话，这里是证明：
本质上是因为 \(f(n)=\mu(n)g(n)\) 且 \(g\) 为完全积性函数。
所以通法是再卷一个 \(g\)，就可以把 \(g\) 全部配成 \(g(n)\) 然后提到卷积外面了。
总之这样就做完了。时间复杂度 \(O(n^{2/3}+k\log k)\)。
2. 杜教筛的 dp 实现 & 整除分块优化
我们先从头推一遍杜教筛的核心式子。
一般是将它写成 \(\displaystyle\sum_{i=1}^nf(i)Sg(n/i)\)，然后用整除分块计算。
这里我们换一种做法。
注意到 \(ij\le n\) 表明 \(i\le \sqrt n\) 或 \(j\le \sqrt n\)。
于是我们可以通过简单的容斥得到：
\(\sqrt n\) 可能不是整数，但 \(i\) 一定是整数，所以取求和上界为 \(B=\lfloor\sqrt n\rfloor\) 即可。
本题中 \(f*g=\varepsilon\)，简单整理一下式子得到：
这个式子可以直接计算，不再需要整除分块。
接下来是对杜教筛本身的优化。
考虑熟知性质 \((a/b)/c=a/(bc)\)，可知我们只需用到 \(Sf(1)\sim Sf(B),Sf(n/B)\sim Sf(n/1)\) 的值。
为了方便，称之为 \(f\) 的块筛。
块筛可以用两个长为 \(B\) 的数组拼起来存储。这个可以看代码实现。
先用线性筛预处理 \(Sf(1)\sim Sf(n^{2/3})\)。
对于块筛的剩余部分，只需应用上面的式子，依次递推求出。（这部分可视作 dp）
过程中还要用到 \(g\) 的块筛。这个显然可以 \(O(\sqrt n)\) 预处理。
时间复杂度为 \(\displaystyle O\left(n^{2/3}+\sum_{i=1}^{n^{1/3}}\sqrt{\dfrac ni}\right)=O(n^{2/3})\)，与原版杜教筛相同。（但常数更优）
进一步地，使用这种写法的杜教筛还可以结合其它科技，使复杂度优化至 \(O\left(\dfrac{n^{2/3}}{\log n}\right)\)。
但这里空白太小，写不下了。有兴趣可以自行了解。
3. 代码实现
除了上面说的优化以外，还对整除运算本身进行了优化。（ dv 函数及其相关部分）
常数很小。目前是最优解 rk1。
下面给出完整代码，仅供参考。
#include<bits/stdc++.h>
#define C const
#define readc(x) (scanf("%d",&x),x)
#define rep(i,l,r) for (int i=l; i<=r; ++i)
typedef long long ll;
typedef __int128 Lll;
using namespace std;
C int N=1e6+8,M=32000;
C int n=readc(n),m=readc(m),K=readc(K);
C int Mx=max(n,m),mx=pow(Mx,2./3);
int k,t; ll s,inv[M]; bitset<M>g; int sg[M],p[80000];
bitset<N>v; char f[N]; int sf[N],Sg[M],Sn[M],Sm[M];
inline int dv(Lll x,int y) { return x*inv[y]>>56; } //返回 x/y 的值
inline int mo(int x,int y) { return x-dv(x,y)*y; }
int gcd(int x,int y) { return y?gcd(y,mo(x,y)):x; }
void prework() {
C int sq=max(K,(int)sqrt(Mx));
C ll base=(1ll<<56);
rep(i,1,sq) inv[i]=base/i+1;
rep(i,1,K) sg[i]=sg[i-1]+(g[i]=(gcd(K,i)==1));
rep(i,K+1,sq) sg[i]=sg[i-1]+(g[i]=g[i-K]);
f[1]=1;
rep(i,2,mx) { //线性筛预处理一部分 f
if (!v[i]) f[i]=(g[mo(i,K)]?-1:0),p[++k]=i;
for (int j=1; j<=k&&(t=i*p[j])<=mx; ++j) {
v[t]=1;
if (!mo(i,p[j])) break;
f[t]=f[i]*f[p[j]];
}
}
rep(i,1,mx) sf[i]=sf[i-1]+f[i];
}
C int phiK=(prework(),sg[K]);
void calc(C int n,int *S) { //用杜教筛求 f 的块筛的第二部分，并存入 S
C int sq=sqrt(n);
rep(i,n/mx+1,sq) S[i]=sf[dv(n,i)]; //已经预处理的部分
rep(i,1,sq) t=dv(n,i),Sg[i]=phiK*dv(t,K)+sg[mo(t,K)]; //求 g 的块筛
for (int i=n/mx; i; --i) { // dp 实现杜教筛
C int n2=dv(n,i),B=sqrt(n2),mid=dv(sq,i);
int s=1+sg[B]*sf[B]-Sg[i];
rep(j,2,mid) s-=g[j]*S[i*j]+f[j]*Sg[i*j];
rep(j,mid+1,B) t=dv(n2,j),s-=g[j]*sf[t]+f[j]*sg[t];
S[i]=s;
}
}
void solve() { //整除分块计算答案
C int Mn=min(n,m),sq=min(n,(int)sqrt(Mx)),Sq=sqrt(m);
rep(i,1,sq) s+=ll(dv(n,i))*f[i]*(i<=Sq?Sg[i]:sg[dv(m,i)]);
if (sq==Mn) { printf("%lld\n",s); return ; }
int t1=dv(n,sq),r1=dv(n,t1);
int t2=dv(m,sq),r2=dv(m,t2);
for (int l=sq+1,lst=sf[sq],cur; l<=Mn; l=min(r1,r2)+1) {
if (l>r1) r1=dv(n,--t1);
if (l>r2) r2=dv(m,--t2);
cur=(r1<r2?Sn[t1]:Sm[t2]);
s+=t1*ll(cur-lst)*sg[t2],lst=cur;
}
printf("%lld\n",s);
}
main() { calc(n,Sn),calc(m,Sm),solve(); }