[atARC123F]Insert Addition

前置知识

下面，先来介绍一下Stern-Brocot Tree的结构：

其是一棵满二叉树，每一个节点都是一个最简分数，其中根为$\frac{1}{1}$

假设前$i$层的中序遍历分数依次为$\frac{y_{i,j}}{x_{i,j}}$（其中$i\ge 1,j\in [1,2^{i})$，即根为第一层），并令$(x_{i,0},y_{i,0})=(1,0)$且$(x_{i,2^{i}},y_{i,2^{i}})=(0,1)$，那么第$i+1$层的从左到右第$j$个分数为$\frac{y_{i,j-1}+y_{i,j}}{x_{i,j-1}+x_{i,j}}$（其中$j\in [1,2^{i}]$）

（特别的，根据定义在$i=0$时，$\{(x_{i},y_{i})\}$即仅包含$(1,0)$和$(0,1)$）

根据定义，不难得到有$f(x_{i})=x_{i+1}$且$f(y_{i})=y_{i+1}$，这可以方便理解下面的结论——

性质1：$\forall i\ge 0,\frac{y_{i,j}}{x_{i,j}}$构成一个严格单调递增的序列（特别的，定义$\frac{1}{0}=\infty$）

性质2：若$x_{i},y_{i}\in N^{*}$，存在$i\ge 0,j\in [0,2^{i})$满足$(x_{i,j},y_{i,j})=(x_{1},y_{1})$且$(x_{i,j+1},y_{i,j+1})=(x_{2},y_{2})$当且仅当$x_{1}y_{2}-x_{2}y_{1}=1$

性质3：若$x,y\in N^{*}$，存在$i\ge 0,j\in [0,2^{i}]$满足$(x_{i,j},y_{i,j})=(x,y)$当且仅当$\gcd(x,y)=1$

性质1和性质2可以归纳得到，性质3可以根据性质2得到，证明过程略

题解

性质4：$\forall i\ge 0,j\in [0,2^{i}]$，有$x_{i,j}+y_{i,j}\ge i+1$

性质5：$\forall i\ge 0,j\in [0,2^{i}]$，有$F^{(i)}(\{a,b\})_{j}=ax_{i,j}+by_{i,j}$

这两个性质同样可以归纳得到，证明过程略

结合性质1和3，不难得到$\{B_{i}\}=\{ax+by\mid x,y\in N^{*},\gcd(x,y)=1$且$ax+by\le N \}$，并且对应的顺序即按照$\frac{y}{x}$从小到大排序，由此即有一个$o(N^{2}\log N)$的暴力做法

考虑优化，不妨假设$\gcd(a,b)=1$（将$N$除以$\gcd(a,b)$即可），那么根据裴蜀定理，存在$ab'-a'b=1$，并根据通解的形式调整使得$a'\in [0,a)$，进而有$b'=\frac{a'b+1}{a}\in [1,b]$

此时，不妨将二元组$(x,y)$变为$(a'x+b'y,ax+by)$，并对条件分析——

性质6：$\gcd(x,y)=\gcd(a'x+b'y,ax+by)$

性质7：$\frac{y_{1}}{x_{1}}<\frac{y_{2}}{x_{2}}$等价于$\frac{a'x_{1}+b'y_{1}}{ax_{1}+by_{1}}<\frac{a'x_{2}+b'y_{2}}{ax_{2}+by_{2}}$

性质6可以归纳并结合辗转相减法得到，性质7直接展开可得，证明过程略

另外，关于$x,y\in N^{*}$的条件，将$x$和$y$用新二元组的$(x,y)$表示，解得即$x,y\in N^{*}$且$\frac{a'}{a}\le \frac{x}{y}\le \frac{b'}{b}$

综上，问题即求分母在$[1,N]$中、值在$[\frac{a'}{a},\frac{b'}{b}]$中的从小到大第$[l,r]$个最简分数（的分母），此时只需要在之前的Stern-Brocot Tree上二分即可

在过程中，即需要查询子树大小，也即求出值在$[\frac{x_{1}}{y_{1}},\frac{x_{2}}{y_{2}}]$中且分母$\in[1,N]$的分数个数

关于这个问题，先将范围差分（即仅考虑$\le \frac{x_{2}}{y_{2}}$），并对$\gcd(x,y)=1$的条件容斥，即令$g(n)$为忽略此条件且分母$\in [1,n]$的分数个数，那么不难得到答案即$\sum_{d=1}^{N}\mu(d)g(\lfloor\frac{N}{d}\rfloor)$

预处理$\mu$的前缀和后，数论分块即变为求$\sqrt{N}$个$g(n)$的值，则有$g(n)=\sum_{y=1}^{n}\left(\lfloor\frac{y\cdot x_{2}}{y_{2}}\rfloor+1\right)$，可以通过类欧几里得算法计算，时间复杂度为$o(\log n)$

注意到每层只有$o(1)$次上述查询，因此总共只有$o(N)$次查询

总复杂度即$o(N\sqrt{N}\log N+len)$（其中$len=r-l+1$），可以通过

（实际上跑不满，常数极小，甚至于一些比较暴力的代码都可以通过）

#include<bits/stdc++.h>
#include<atcoder/math>
using namespace std;
#define N 300005 
#define ll long long
#define pii pair<int,int>
#define fi first
#define se second
vector<int>ans;
int n,aa,bb,vis[N],p[N],mu[N];
ll l,r;
pii a,b;
bool cmp(pii x,pii y){
    return (ll)x.fi*y.se<(ll)x.se*y.fi;
}
ll get_g(pii k,int n){
    return atcoder::floor_sum(n,k.se,k.fi,k.fi+k.se);
}
ll get_tot(pii k){
    if (cmp(k,a))k=a;
    if (cmp(b,k))k=b;
    ll ans=0;
    for(int i=1,j;i<=n;i=j+1){
        j=n/(n/i);
        ans+=(ll)(mu[j]-mu[i-1])*get_g(k,n/i);
    }
    return ans;
}
void dfs(pii x,pii y){
    if ((cmp(b,x))||(cmp(y,a)))return;
    pii m=make_pair(x.fi+y.fi,x.se+y.se);
    if (m.se>n)return;
    dfs(x,m);
    if ((cmp(a,m))&&(cmp(m,b)))ans.push_back(m.se);
    dfs(m,y);
}
void calc(pii x,pii y,ll l,ll r,ll s){
    if ((cmp(b,x))||(cmp(y,a))||(l>s)||(r<1))return;
    if ((l<=1)&&(s<=r)){
        dfs(x,y);
        return;
    }
    pii m=make_pair(x.fi+y.fi,x.se+y.se);
    ll ss=get_tot(m)-get_tot(x);
    calc(x,m,l,r,ss-1);
    if ((l<=ss)&&(ss<=r)&&(cmp(a,m))&&(cmp(m,b)))ans.push_back(m.se);
    calc(m,y,l-ss,r-ss,s-ss);
}
int main(){
    mu[1]=1;
    for(int i=2;i<N;i++){
        if (!vis[i]){
            p[++p[0]]=i;
            mu[i]=-1;
        }
        for(int j=1;(j<=p[0])&&(i*p[j]<N);j++){
            vis[i*p[j]]=1;
            if (i%p[j])mu[i*p[j]]=mu[i]*mu[p[j]];
            else{
                mu[i*p[j]]=0;
                break;
            }
        }
    }
    for(int i=1;i<N;i++)mu[i]+=mu[i-1];
    scanf("%d%d%d%lld%lld",&aa,&bb,&n,&l,&r);
    int g=__gcd(aa,bb);
    aa/=g,bb/=g,n/=g;
    for(int i=0;i<aa;i++)
        if (((ll)i*bb+1)%aa==0){
            a=make_pair(i,aa);
            b=make_pair(((ll)i*bb+1)/aa,bb);
            break;
        }
    if (l>1)l--;
    else ans.push_back(a.se);
    r--;
    pii x=make_pair(0,1),y=make_pair(1,0);
    calc(x,y,l,r,get_tot(y)-get_tot(x)-1);
    if (get_tot(y)-get_tot(x)==r)ans.push_back(b.se);
    for(int i=0;i<ans.size();i++)ans[i]*=g;
    printf("%d",ans[0]);
    for(int i=1;i<ans.size();i++)printf(" %d",ans[i]);
    printf("\n");
    return 0;
}