@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 的既约分数可以用如下式子计算：
\[ans = \sum_{i=1}^{n}\sum_{j=1}^{\lfloor x*i\rfloor}[gcd(i, j) = 1]\]

对其进行反演可得：
\[ans = \sum_{i=1}^{n}\sum_{d=1}^{d|i}\mu(d)*\lfloor\frac{x*i}{d}\rfloor \\ = \sum_{a*b\le n}\mu(a)*\lfloor x*b \rfloor \\ = \sum_{i=1}^{n}\mu(i)*\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\lfloor x*j\rfloor\]

后面那一项如果以 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)。

@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@

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

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

转载于:https://www.cnblogs.com/Tiw-Air-OAO/p/11236142.html