【洛谷1587】[NOI2016] 循环之美（杜教筛）

点此看题面

• 所有\(x\in[1,n],y\in[1,m]\)中，共有多少个值不相同的分数\(\frac xy\)，满足在\(k\)进制下小数点后是纯循环的（包括整数）。
• \(n,m\le10^9,k\le2000\)

计算式的表示

首先要求值不相同，等价于只考虑最简分数，也就是说一个必要条件是\(gcd(i,j)=1\)。

然后考虑什么样的数小数点后是纯循环的，根据我们的生活经验，当\(k=10\)时的充要条件为\(j\)不是\(2\)或\(5\)的倍数，因此可以猜想一般情况下的条件就是\(gcd(j,k)=1\)。

严谨的证明是设循环节长度为\(l\)，则：

\[\frac ij-\lfloor\frac ij\rfloor=\frac ij\times k^l-\lfloor\frac ij\times k^l\rfloor\\ i-\lfloor\frac ij\rfloor\times j=i\times k^l-\lfloor\frac ij\times k^l\rfloor\times j\\ i\equiv i\times k^l(mod\ j) \]

因为\(gcd(i,j)=1\)，可以进一步化简得到：

\[k^l\equiv1(mod\ j) \]

由于存在这样一个\(l\)，因此得到：

\[gcd(j,k)=1 \]

所以我们列出此题的计算式：

\[\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1] \]

递归求解+杜教筛处理小情况

考虑把\([gcd(j,k)=1]\)一项拆开，得到：

\[\sum_{d|k}\mu(d)\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac md\rfloor}[gcd(i,jd)=1] \]

由于\(gcd(i,jd)=1\Leftrightarrow gcd(i,j)=1\wedge gcd(i,d)=1\)，所以可以进一步写成：

\[\sum_{d|k}\mu(d)\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac md\rfloor}[gcd(i,j)=1][gcd(i,d)=1] \]

观察发现后面同样有两个\(gcd\)，且两个\(gcd\)中同样有一元相同。

因此我们设\(f(n,m,k)=\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1]\)，由此得到：

\[f(n,m,k)=\sum_{d|k}\mu(d) f(\lfloor\frac md\rfloor,n,d) \]

边界条件是\(n=0\)或\(m=0\)时\(f(n,m,k)=0\)。

注意到若枚举到的\(d=1\)则\(n,m\)只是简单交换而不会变小。

不过此时问题就退化成求\(\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=1]\)。

这是一个经典的莫比乌斯反演问题，直接给出转化后的式子为\(\sum_{i=1}^{\min\{n,m\}}\mu(i)\lfloor\frac ni\rfloor\lfloor\frac mi\rfloor\)。

由于\(n,m\)范围很大，需要杜教筛。

代码：\(O(\sqrt nlog^2n+n^{\frac23})\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define S 1000000
#define LL long long
using namespace std;
int n,m,k;
namespace Du//杜教筛
{
	int Pt,P[S+5],mu[S+5],s[S+5];I void Sieve()//线性筛预处理小范围的μ
	{
		mu[1]=1;for(RI i=2,j;i<=S;++i) for(!P[i]&&(mu[P[++Pt]=i]=-1),
			j=1;j<=Pt&&i*P[j]<=S;++j) if(P[i*P[j]]=1,i%P[j]) mu[i*P[j]]=-mu[i];else break;
		for(RI i=1;i<=S;++i) s[i]=s[i-1]+mu[i];
	}
	map<int,int> p;I int SMu(CI x)//杜教筛
	{
		if(x<=S) return s[x];if(p.count(x)) return p[x];//小范围直接调用；记忆化
		RI t=1;for(RI l=2,r;l<=x;l=r+1) r=x/(x/l),t-=(r-l+1)*SMu(x/l);return p[x]=t;//除法分块
	}
}
struct Data {int n,m,k;I bool operator < (Con Data& o) Con {return n^o.n?n<o.n:(m^o.m?m<o.m:k<o.k);}};
map<Data,LL> ans;I LL F(CI n,CI m,CI k)//递归求解
{
	if(!n||!m) return 0;if(ans.count((Data){n,m,k})) return ans[(Data){n,m,k}];
	LL t=0;if(k==1)//如果k=1
	{
		for(RI l=1,r;l<=min(n,m);l=r+1)//除法分块
			r=min(n/(n/l),m/(m/l)),t+=1LL*(n/l)*(m/l)*(Du::SMu(r)-Du::SMu(l-1));
		return ans[(Data){n,m,k}]=ans[(Data){m,n,k}]=t;
	}
	#define Calc(x) (Du::mu[x]?Du::mu[x]*F(m/(x),n,x):0)//一个优化，μ(x)=0时无需递归
	for(RI i=1;i*i<=k;++i) !(k%i)&&(t+=Calc(i),i^(k/i)&&(t+=Calc(k/i)));return ans[(Data){n,m,k}]=t;//枚举因数
}
int main()
{
	return Du::Sieve(),scanf("%d%d%d",&n,&m,&k),printf("%lld\n",F(n,m,k)),0;
}