[NOI2016] 循环之美

description

求

[sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=1][frac ij在k进制下是纯循环小数] ]

data range

(n,mle 10^9,kle 2000)

solution

(frac ij)在(k)进制下为纯循环小数当且仅当(gcd(j,k)=1)

证明戳这里

那么原式

[=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1]\ =sum_{j=1}^m[gcd(j,k)=1]sum_{i=1}^n[gcd(i,j)=1]\ =sum_{j=1}^m[gcd(j,k)=1]sum_{i=1}^nsum_{dmid i,dmid j}mu(d)\ =sum_{d=1}^nmu(d)sum_{dmid j}^m[gcd(j,k)=1]sum_{dmid i}^n1\ =sum_{d=1}^nmu(d)sum_{j=1}^{lfloor frac md floor}[gcd(jd,k)=1]lfloor frac nd floor\ =sum_{d=1}^nlfloor frac nd floormu(d)[gcd(d,k)=1]sum_{j=1}^{lfloor frac md floor}[gcd(j,k)=1] ]

不妨设

[F(n)=sum_{i=1}^n[gcd(i,k)=1]\ S(n,k)=sum_{i=1}^nmu(i)[gcd(i,k)=1]\ ]

那么原式

[=sum_{d=1}^nlfloor frac nd floor F(lfloor frac md floor)(S(d,k)-S(d-1,k)) ]

显然如果可以快速求出(S(n,k),F(n)) ，那么通过整除分块就可以快速求得原式了

先来康康(F(n)),我们有

[F(n)=lfloorfrac nk floor F(k)+F(nmod k) ]

只用预处理处(nle k)的(F(n))就可以做到快速查询了

再来康康(S(n,k))

[S(n,k)=sum_{i=1}^nmu(i)[gcd(i,k)=1]\ =sum_{i=1}^nmu(i)sum_{dmid i,dmid k}mu(d)\ =sum_{dmid k}mu(d)sum_{dmid i}^nmu(i)\ =sum_{dmid k}mu(d)sum_{i=1}^{lfloor frac nd floor}mu(id) ]

显然若(gcd(i,d)>1) 则(mu(id)=0)

因此

[=sum_{dmid k}mu(d)sum_{i=1}^{lfloor frac nd floor}mu(id)[gcd(i,d)=1]\ =sum_{dmid k}mu(d)sum_{i=1}^{lfloor frac nd floor}mu(i)mu(d)[gcd(i,d)=1]\ =sum_{dmid k}mu(d)^2sum_{i=1}^{lfloor frac nd floor}mu(i)[gcd(i,d)=1]\ =sum_{dmid k}mu(d)^2S(lfloor frac nd floor,d) ]

递归求解即可（记得记忆化）

注意递归边界

(n=0) 时,(S(0,k)=0)

(k=1) 时,(S(n,1)=sum_{i=1}^nmu(i)[gcd(1,k)=1]=sum_{i=1}^nmu(i))

杜教筛之即可

完结撒花。。。

time complexity

(mathcal O(能过))

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+5,Bas=2001;
bool flag[N];int pr[N],pcnt,mu[N],smu[N],f[N];
unordered_map<ll,ll>mp;
inline ll qs(int a,int b){return 1ll*a*Bas+b;}
ll ss(int n,int k)
{
	if(!n)return 0;
	ll now=qs(n,k); 
	if(mp.count(now))return mp[now];
	if(k==1)
	{
		if(n<=1e6)return smu[n];
		ll ret=1;
		for(int l=2,r;l<=n;l=r+1)
		{
			r=n/(n/l);
			ret-=ss(n/l,1)*(r-l+1);
		}return mp[now]=ret;
	}ll ret=0;
	for(int i=1;i*i<=k;++i)
	{
		if(k%i)continue;
		if(mu[i])ret+=ss(n/i,i);
		if(i*i!=k&&mu[k/i])ret+=ss(n/(k/i),k/i);
	}return mp[now]=ret;
}
int k;
int gcd(int x,int y){return y?gcd(y,x%y):x;}
inline void pre(int n)
{
	mu[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!flag[i])pr[++pcnt]=i,mu[i]=-1;
		for(int j=1;j<=pcnt;++j)
		{
			int num=i*pr[j];if(num>n)break;
			flag[num]=1;
			if(i%pr[j])mu[num]=-mu[i];
			else{mu[num]=0;break;}
		}
	}
	for(int i=1;i<=n;++i)smu[i]=smu[i-1]+mu[i];
	for(int i=1;i<=k;++i)f[i]=f[i-1]+(gcd(i,k)==1);
}
inline int sf(int n){return n/k*f[k]+f[n%k];}
int n,m;
int main()
{
	scanf("%d%d%d",&n,&m,&k);pre(1e6);
	int mn=min(n,m);ll ans=0;
	for(int l=1,r;l<=mn;l=r+1)
	{
		r=min(n/(n/l),m/(m/l));
		ans+=(ss(r,k)-ss(l-1,k))*(n/l)*sf(m/l);
	}printf("%lld
",ans);
	return 0;
}