【洛谷1587】[NOI2016] 循环之美(杜教筛)

点此看题面

  • 所有(xin[1,n],yin[1,m])中,共有多少个值不相同的分数(frac xy),满足在(k)进制下小数点后是纯循环的(包括整数)。
  • (n,mle10^9,kle2000)

计算式的表示

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

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

严谨的证明是设循环节长度为(l),则:

[frac ij-lfloorfrac ij floor=frac ij imes k^l-lfloorfrac ij imes k^l floor\ i-lfloorfrac ij floor imes j=i imes k^l-lfloorfrac ij imes k^l floor imes j\ iequiv i imes k^l(mod j) ]

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

[k^lequiv1(mod j) ]

由于存在这样一个(l),因此得到:

[gcd(j,k)=1 ]

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

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

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

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

[sum_{d|k}mu(d)sum_{i=1}^nsum_{j=1}^{lfloorfrac md floor}[gcd(i,jd)=1] ]

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

[sum_{d|k}mu(d)sum_{i=1}^nsum_{j=1}^{lfloorfrac md floor}[gcd(i,j)=1][gcd(i,d)=1] ]

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

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

[f(n,m,k)=sum_{d|k}mu(d) f(lfloorfrac md floor,n,d) ]

边界条件是(n=0)(m=0)(f(n,m,k)=0)

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

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

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

由于(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
",F(n,m,k)),0;
}
败得义无反顾,弱得一无是处
原文地址:https://www.cnblogs.com/chenxiaoran666/p/Luogu1587.html