[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$),可以通过
(实际上跑不满,常数极小,甚至于一些比较暴力的代码都可以通过)
1 #include<bits/stdc++.h> 2 #include<atcoder/math> 3 using namespace std; 4 #define N 300005 5 #define ll long long 6 #define pii pair<int,int> 7 #define fi first 8 #define se second 9 vector<int>ans; 10 int n,aa,bb,vis[N],p[N],mu[N]; 11 ll l,r; 12 pii a,b; 13 bool cmp(pii x,pii y){ 14 return (ll)x.fi*y.se<(ll)x.se*y.fi; 15 } 16 ll get_g(pii k,int n){ 17 return atcoder::floor_sum(n,k.se,k.fi,k.fi+k.se); 18 } 19 ll get_tot(pii k){ 20 if (cmp(k,a))k=a; 21 if (cmp(b,k))k=b; 22 ll ans=0; 23 for(int i=1,j;i<=n;i=j+1){ 24 j=n/(n/i); 25 ans+=(ll)(mu[j]-mu[i-1])*get_g(k,n/i); 26 } 27 return ans; 28 } 29 void dfs(pii x,pii y){ 30 if ((cmp(b,x))||(cmp(y,a)))return; 31 pii m=make_pair(x.fi+y.fi,x.se+y.se); 32 if (m.se>n)return; 33 dfs(x,m); 34 if ((cmp(a,m))&&(cmp(m,b)))ans.push_back(m.se); 35 dfs(m,y); 36 } 37 void calc(pii x,pii y,ll l,ll r,ll s){ 38 if ((cmp(b,x))||(cmp(y,a))||(l>s)||(r<1))return; 39 if ((l<=1)&&(s<=r)){ 40 dfs(x,y); 41 return; 42 } 43 pii m=make_pair(x.fi+y.fi,x.se+y.se); 44 ll ss=get_tot(m)-get_tot(x); 45 calc(x,m,l,r,ss-1); 46 if ((l<=ss)&&(ss<=r)&&(cmp(a,m))&&(cmp(m,b)))ans.push_back(m.se); 47 calc(m,y,l-ss,r-ss,s-ss); 48 } 49 int main(){ 50 mu[1]=1; 51 for(int i=2;i<N;i++){ 52 if (!vis[i]){ 53 p[++p[0]]=i; 54 mu[i]=-1; 55 } 56 for(int j=1;(j<=p[0])&&(i*p[j]<N);j++){ 57 vis[i*p[j]]=1; 58 if (i%p[j])mu[i*p[j]]=mu[i]*mu[p[j]]; 59 else{ 60 mu[i*p[j]]=0; 61 break; 62 } 63 } 64 } 65 for(int i=1;i<N;i++)mu[i]+=mu[i-1]; 66 scanf("%d%d%d%lld%lld",&aa,&bb,&n,&l,&r); 67 int g=__gcd(aa,bb); 68 aa/=g,bb/=g,n/=g; 69 for(int i=0;i<aa;i++) 70 if (((ll)i*bb+1)%aa==0){ 71 a=make_pair(i,aa); 72 b=make_pair(((ll)i*bb+1)/aa,bb); 73 break; 74 } 75 if (l>1)l--; 76 else ans.push_back(a.se); 77 r--; 78 pii x=make_pair(0,1),y=make_pair(1,0); 79 calc(x,y,l,r,get_tot(y)-get_tot(x)-1); 80 if (get_tot(y)-get_tot(x)==r)ans.push_back(b.se); 81 for(int i=0;i<ans.size();i++)ans[i]*=g; 82 printf("%d",ans[0]); 83 for(int i=1;i<ans.size();i++)printf(" %d",ans[i]); 84 printf("\n"); 85 return 0; 86 }