@hdu

目录

• @description@
• @solution@
• @accepted code@
• @details@

---

@description@

询问第 k 小的分子分母 ≤ n 的既约分数。

Input
第一行包含一个整数 T(T≤10^2)，表示数据组数。
接下来 T 行每行两个整数 n, k，表示一组询问。保证询问的答案在 (0,1] 范围内。

Output
输出 T 行，每行包含一个分数 p/q，要求 gcd(p, q) = 1。

Sample Input
5
4 6
5 1
9 9
3 4
7 11
Sample Output
1/1
1/5
1/3
1/1
3/5

Hint
杭电支持 __int128。

@solution@

寻找第 k 小，不难想到使用二分法。
则假如二分到一个值 x，≤ x 的既约分数可以用如下式子计算：

[ans = sum_{i=1}^{n}sum_{j=1}^{lfloor x*i floor}[gcd(i, j) = 1] ]

对其进行反演可得：

[ans = sum_{i=1}^{n}sum_{d=1}^{d|i}mu(d)*lfloorfrac{x*i}{d} floor \ = sum_{a*ble n}mu(a)*lfloor x*b floor \ = sum_{i=1}^{n}mu(i)*sum_{j=1}^{lfloorfrac{n}{i} floor}lfloor x*j floor]

后面那一项如果以 j 为横坐标，x 为常数，可以看成直线 y = x*j 下的整点数量，联想到类欧几里得算法。
但是 x 是个实数的话类欧几里得是不可行的。不过假如我们二分直接使用分数进行二分的话，就的确可以使用类欧几里得算法。
加上预处理莫比乌斯函数前缀和 + 整除分块的话，可以在 (O(sqrt{n}*log n)) 的时间内完成一次计算。

实数二分需要一个迭代次数，那么我们多少次才能保证精度呢？
假如 a/b 与 c/d 都是两个分子分母 ≤ n 的不相等的分数，则 a/b - c/d = (ad - bc)/(b*d) 最小时，分子 ad - bc = 1，分母 b*d = n*n。
所以精度的一个上限是 1/(n*n)，于是取 log(n*n) 次二分就可以保证精度。这样最大次数为 40。

下一个问题是，我们如何找到 > p/q 且最小的既约分数呢？
我们可以在 SB tree（Stern-Brocot Tree）上进行寻找。这棵树在具体数学上有介绍它的性质，WC2019 上也有讲一点。
这里有一篇简略的博客。

我们这道题只需要知道 SB Tree 的生成方式，以及它的几个性质：
（1）SB Tree 是一个既约分数的二叉搜索树。
（2）SB Tree 父亲的分母小于儿子的分母。
而这些足以让我们实现我们的目的。寻找答案的这个过程是 O(n)。

update in 2020/07/05：我试了一下，貌似不用SB-Tree也可以找。枚举分母然后算最小合法分子即可，一样是线性。不过由于需要做除法所以要慢一些。

@accepted code@

#pragma comment(linker, "/STACK:102400000,102400000")
#include<cstdio>
#include<cmath>
typedef long long ll;
const int MAXN = 1000000;
__int128 gcd(__int128 a, __int128 b) {return (b == 0) ? a : gcd(b, a % b);}
struct frac{
	__int128 a, b;
	frac(__int128 _a=0, __int128 _b=0):a(_a), b(_b) {}
	friend bool operator < (frac a, frac b) {return a.a*b.b < b.a*a.b;}
	friend frac operator + (frac a, frac b) {
		__int128 d = gcd(a.b*b.a + a.a*b.b, a.b*b.b);
		return frac((a.b*b.a + a.a*b.b)/d, a.b*b.b/d);
	}
	friend frac operator / (frac a, __int128 k) {return frac(a.a, k*a.b);}
};
int n; ll k;
bool search_answer(frac p) {
	frac ans = frac(1, 1), a = frac(0, 1), b = frac(1, 1);
	while( a.b + b.b <= n ) {
		frac c = frac(a.a + b.a, a.b + b.b);
		if( p < c )
			ans = b = c;
		else a = c;
	}
	printf("%d/%d
", int(ans.a), int(ans.b));
}
ll smiu[MAXN + 5];
__int128 func(__int128 A, __int128 B, __int128 C, __int128 n) {
	if( A/C ) return func(A % C, B, C, n) + A/C*(n + 1)*n/2;
	if( B/C ) return func(A, B % C, C, n) + B/C*(n + 1);
	if( A == 0 ) return 0;
	__int128 m = (A*n + B) / C;
	if( m == 0 ) return 0;
	return m*n - func(C, C - 1 - B, A, m - 1);
}
ll check(frac x) {
	ll ans = 0; int p = 1, q = 1;
//	printf("? %d %d
", int(x.a), int(x.b));
	while( p <= n ) {
		q = p, p = (n/(n/q));
		ll s = smiu[p] - smiu[q - 1];
		ans += s*func(x.a, 0, x.b, n/p);
	//	printf("! %lld
", ans);
		p++;
	}
	return ans;
}
void solve() {
	scanf("%d%lld", &n, &k);
	int lim = log2(1LL*n*n) + 1;
	frac le = frac(0, 1), ri = frac(1, 1);
	for(int i=0;i<=lim;i++) {
		frac mid = (le + ri) / 2;
		if( check(mid) < k )
			le = mid;
		else ri = mid;
	}
	search_answer(le);
}
bool nprm[MAXN + 5];
void init() {
	for(int i=1;i<=MAXN;i++)
		smiu[i] = 1;
	for(int i=2;i<=MAXN;i++) {
		if( !nprm[i] ) {
			for(int j=i;j<=MAXN;j+=i) {
				nprm[j] = true;
				if( (j/i) % i == 0 )
					smiu[j] = 0;
				else smiu[j] *= -1;
			}
		}
	}
	for(int i=1;i<=MAXN;i++)
		smiu[i] += smiu[i-1];
}
int main() {
	init();
	int T; scanf("%d", &T);
	while( T-- ) solve();
}

@details@

算是还不错的一道题，不过涉及到的知识点比较冷门。。。
我才不会说理解类欧几里得算法花了我一个晚上的时间。

找答案的时候还不能在树上递归找，不然会爆栈。正确的操作是迭代找答案。