@hdu - 6584@ Meteor
@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 的既约分数可以用如下式子计算：
对其进行反演可得：
后面那一项如果以 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\n", 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\n", 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\n", 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@
算是还不错的一道题，不过涉及到的知识点比较冷门。。。
我才不会说理解类欧几里得算法花了我一个晚上的时间。
找答案的时候还不能在树上递归找，不然会爆栈。正确的操作是迭代找答案。
