P8058 [BalkanOI2003] Farey 序列 题解

分析

考虑二分答案。

对于当前二分的答案 \(x\)，设 \(cnt\) 表示 Farey 序列中 \(\frac{p}{q} \le x\) 的满足条件的数量。对于一组 \((i,j)\)，若 \(\frac{j}{i}\le x\)，则 \(j \le\lfloor i \times x \rfloor\)。得到暴力式子：

\[cnt=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{\lfloor i \times x \rfloor}[\gcd(i,j)=1] \]

定义 \(f_i=\sum\limits_{j=1}^{\lfloor i \times x \rfloor}[\gcd(i,j)=1]\)。根据容斥，有：

\[f_i=\lfloor i \times x \rfloor-\sum\limits_{j=1}^{\lfloor i \times x \rfloor}[\gcd(i,j) >1] \]

将 \(\gcd(i,j)\) 的值表示出来：

\[f_i=\lfloor i \times x \rfloor-\sum\limits_{j=1}^{\lfloor i \times x \rfloor}[\gcd(i,j) =k \land k>1] \]

把 \(k\) 提出来：

\[f_i=\lfloor i \times x \rfloor-\sum\limits_{k=2}^{i}\sum\limits_{j=1}^{\lfloor\lfloor\frac{ i }{k}\rfloor \times x\rfloor}[\gcd(i,j) =1] \]

不难发现，第二个求和等价于 \(f_{\lfloor \frac{i}{k} \rfloor}\)，则有：

\[f_i=\lfloor i \times x \rfloor-\sum\limits_{k=1}^{i}f_{\lfloor \frac{i}{k} \rfloor} \]

化简得：

\[f_i=\lfloor i \times x \rfloor-\sum\limits_{k|i \land 1<k<i}^{}f_k \]

\(cnt=\sum\limits_{i=1}^{n}f_i\)。枚举 \(i\) 的约数，则能够在 \(O(n\sqrt{n})\) 的复杂度求出 \(cnt\)，然后复杂度就是 \(O(n \sqrt{n} \log n)\)。

对于求 \(p,q\) 的值，暴力枚举即可。

代码

#include<bits/stdc++.h>
using namespace std;
//#define int long long
#define re register
#define il inline
#define pii pair<int,int>
#define x first
#define y second
#define gc getchar()
#define rd read()
#define debug() puts("------------")

namespace yzqwq{
	il int read(){
		int x=0,f=1;char ch=gc;
		while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=gc;}
		while(ch>='0'&&ch<='9') x=(x<<1)+(x<<3)+(ch^48),ch=gc;
		return x*f;
	}
	il int qmi(int a,int b,int p){
		int ans=1;
		while(b){
			if(b&1) ans=ans*a%p;
			a=a*a%p,b>>=1;
		}
		return ans;
	}
	il auto max(auto a,auto b){return (a>b?a:b);}
	il auto min(auto a,auto b){return (a<b?a:b);}
	il int gcd(int a,int b){
		if(!b) return a;
		return gcd(b,a%b);
	}
	il int lcm(int a,int b){
		return a/gcd(a,b)*b;
	}
	il void exgcd(int a,int b,int &x,int &y){
		if(!b) return x=1,y=0,void(0);
		exgcd(b,a%b,x,y);
		int t=x;
		x=y,y=t-a/b*x;
		return ;
	}
	mt19937 rnd(time(0));
}
using namespace yzqwq;

const int N=4e4+10;
double eps=1e-10;
int n,k,f[N];

il bool check(double x){
	int cnt=0;
	for(re int i=1;i<=n;++i){
		f[i]=i*x;
		for(re int j=2;j*j<=i;++j){
			if(i%j==0){
				f[i]-=f[j];
				if(i/j!=j) f[i]-=f[i/j];
			}
		}
		cnt+=f[i];
	}
	return cnt>=k;
}

il void solve(){
	n=rd,k=rd;
	double l=0.0,r=1.0,ans=0.0;
	while((r-l)>eps){
		double mid=(l+r)/2.0;
		if(check(mid)) ans=mid,r=mid;
		else l=mid;
	}
	int p=0,q=0;
	double c=40001.0;
	for(re int Q=1;Q<=40001;++Q){
		int P=Q*ans;
		if(fabs(double(P*1.0/Q)-ans)<c) c=fabs(double(P*1.0/Q)-ans),p=P,q=Q;
	}
	printf("%lld %lld\n",p,q);
	return ;
}

signed main(){
//	freopen(".in","r",stdin);
//	freopen(".out","w",stdout);
	int t=1;while(t--)
	solve();
	return 0;
}