连分数，Stern-Brocot 树和 Farey 序列 学习笔记

· · 个人记录

连分数，Stern-Brocot 树和 Farey 序列是三种维护有理分数的结构，三者有很强的联系。

连分数

对于数列 a，连分数 [a_0,a_1,\cdots,a_n] 表示展开式

x=a_0+\frac 1{a_1+\frac 1{a_2+\frac 1{\cdots\frac 1{a_n}}}}

此处我们考虑系数全为整数，且 a_{1\sim n}>0 的连分数，这被称为简单连分数（下文简称为连分数）。成 x_k=[a_0,a_1,\cdots,a_k] 为 x 的 k 阶渐进分数，r_k=[a_k,a_{k+1},\cdots,a_n] 为 x 的第 k 个余项。

不难发现，因为 a_{1\sim n}>0，则 r_{1\sim k}>0，所以每个有理数 x=\frac pq 有且仅有两种连分数表示：令 r_0=x,a_k=\lfloor r_k\rfloor,r_{k+1}=\frac 1{r_k-a_k}，当 r 为整数时停止。另一种表示是将最后一项拆成 [a_n-1,1]，前者被称为标准表示。这个过程是辗转相除，因此只会有 \log 项。

\begin{aligned}\frac{11}{45}&=0+\frac{11}{45}\\\frac{45}{11}&=4+\frac 1{11}\\\frac{11}1&=11\\\frac{11}{45}&=[0,4,11]\end{aligned}

考虑一个连分数的渐进分数怎么求。注意到当 a_{0\sim k-1} 确定，x_k=\frac{p_k}{q_k} 显然可以被表示为关于 a_k 的一次函数。设 p_k=A_ka_k+B_k,q_k=C_ka_k+D_k，则有：

\begin{aligned}\frac{p_{k+1}}{q_{k+1}}&=[a_0,a_1,\cdots,a_k,a_{k+1}]\\&=[a_0,a_1,\cdots,a_k+\frac 1{a_{k+1}}]\\&=\frac{A_k(a_k+\frac 1{a_{k+1}})+B_k}{C_k(a_k+\frac 1{a_{k+1}})+D_k}\\&=\frac{(A_ka_k+B_k)A_{k+1}+A_k}{(C_ka_k+D_k)A_{k+1}+C_k}\\\end{aligned}

则 B_{k+1}=A_k,A_{k+1}=p_k。因此 p_k=A_ka_k+B_k=p_{k-1}a_k+A_{k-1}=p_{k-1}a_k+p_{k-2}。同理 q_k=q_{k-1}a_k+q_{k-2}。 反推得边界为 x_{-1}=\frac 10,x_{-2}=\frac 01。

P10788

只考虑 \leq\frac12 的数，\geq2 的数对应这些数的倒数。考虑这些数的生成方式，我们规定禁用 -2，则操作是若 \frac ab\leq\frac12 就加上 2k，否则取倒数。可以证明禁用 -2 无影响，取出第一个被 -2 生成的数，它的上一步是 \frac{a+2b}b\geq2 且被取倒数生成，矛盾。

发现这个过程就是连分数 [a_0,a_1,\cdots,a_n]，要求 a 全为偶数且 a_1\sim a_n 为正。不难证明这个过程生成每个分数的路径唯一且生成的分数最简。直接搜可以 O(ans)。

#include<bits/stdc++.h>
using namespace std;
int n,m;
long long dfs(int a,int b){
  long long ans=(a<=m&&b<=n)+(a<=n&&b<=m);
  for(a+=2*b;a<=max(n,m);a+=2*b)ans+=dfs(b,a);
  return ans;
}
int main(){
  return cin>>n>>m,cout<<dfs(0,1)-2<<'\n',0;
}

考虑优化成一次计算多个连分数。我们给连分数的某个位置留空，钦定最大值（多个最大取最前/后）不填。此时我们可以把搜出来的连分数表示为 \frac{ax+b}{cx+d} 以及 x\geq lim，除一下就知道 x 的个数。剪掉加入最大值 x 后不合法的，可以通过。

#include<bits/stdc++.h>
using namespace std;
int n,m;
long long dfs(int a,int b,int c,int d,int lim,bool type){
  if(1ll*a*lim+b>n)return 0;
  long long ans=type?min(a?(n-b)/a:m,(m-d)/c)-lim+1+max((n-d)/c-lim+1,0):dfs(0,d,2*d,b,lim,1);
  for(int i=1;(type?1ll*(a+2*i*c)*max(lim,i+1)+b+2*i*d:2ll*(b+2*i*d)*max(lim,i)+d)<=m;i++)ans+=dfs(c,d,a+2*i*c,b+2*i*d,max(lim,i+type),type);
  return ans;
}
int main(){
  cin>>n>>m;
  if(n>m)swap(n,m);
  return cout<<dfs(0,0,0,1,1,0)<<'\n',0;
}

P7739

显然 f 就是连分数。考虑用线性变换刻画连分数，设当前分数为 \begin{bmatrix}p\\q\end{bmatrix}，从下往上展开，变为 \frac 1{a_n+\frac pq}=\frac q{a_nq+p}，相当于左乘矩阵 \begin{bmatrix}0&1\\1&a_n\end{bmatrix}。方便起见，a_0 也这么做，最后得到的是倒数。

但是这题维护的是操作序列而不是连分数。考虑将操作序列也看作线性变换。对操作一，从 \begin{bmatrix}0&1\\1&a_n\end{bmatrix} 变为 \begin{bmatrix}0&1\\1&a_n+1\end{bmatrix}，相当于矩阵 \begin{bmatrix}0&1\\1&a_n\end{bmatrix}^{-1}\begin{bmatrix}0&1\\1&a_n+1\end{bmatrix}=\begin{bmatrix}1&1\\0&1\end{bmatrix}；对操作二，当 a_n>1，相当于 \begin{bmatrix}0&1\\1&a_n\end{bmatrix}^{-1}\begin{bmatrix}0&1\\1&a_n-1\end{bmatrix}\begin{bmatrix}0&1\\1&1\end{bmatrix}^2=\begin{bmatrix}0&-1\\1&2\end{bmatrix}；当 a_n=1，相当于 \begin{bmatrix}0&1\\1&1\end{bmatrix}^{-1}\begin{bmatrix}0&1\\1&a_{n-1}\end{bmatrix}^{-1}\begin{bmatrix}0&1\\1&a_{n-1}+1\end{bmatrix}\begin{bmatrix}0&1\\1&1\end{bmatrix}=\begin{bmatrix}0&-1\\1&2\end{bmatrix}，两种情况恰好对应相同的变换。

至此问题变为维护矩阵的乘积。平衡树维护，有反转和翻转的标记，对应的要维护四个信息对应标记的四种状态，下传时交换这些矩阵。

#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int n,m,rt;
string opt,s;
struct matrix{
  int w,x,y,z;
  matrix(){
    w=z=1,x=y=0;
  }
  matrix(int a,int b,int c,int d){
    w=a,x=b,y=c,z=d;
  }
  matrix operator*(matrix a){
    return matrix((1ll*w*a.w+1ll*x*a.y)%mod,(1ll*w*a.x+1ll*x*a.z)%mod,(1ll*y*a.w+1ll*z*a.y)%mod,(1ll*y*a.x+1ll*z*a.z)%mod);
  }
};
const matrix W(1,1,0,1),E(0,mod-1,1,2);
struct node{
  int tr[2],sz,r;
  bool tag1,tag2;
  matrix v0,v1,p00,p01,p10,p11;
};
template<int maxn>struct FHQtreap{
  node tr[maxn];
  int cnt;
  void pushdown1(int pos){
    if(pos)tr[pos].tag1^=1,swap(tr[pos].tr[0],tr[pos].tr[1]),swap(tr[pos].p00,tr[pos].p10),swap(tr[pos].p01,tr[pos].p11);
  }
  void pushdown2(int pos){
    if(pos)tr[pos].tag2^=1,swap(tr[pos].v0,tr[pos].v1),swap(tr[pos].p00,tr[pos].p01),swap(tr[pos].p10,tr[pos].p11);
  }
  void pushdown(int pos){
    if(tr[pos].tag1)pushdown1(tr[pos].tr[0]),pushdown1(tr[pos].tr[1]),tr[pos].tag1=0;
    if(tr[pos].tag2)pushdown2(tr[pos].tr[0]),pushdown2(tr[pos].tr[1]),tr[pos].tag2=0;
  }
  void pushup(int pos){
    tr[pos].sz=tr[tr[pos].tr[0]].sz+tr[tr[pos].tr[1]].sz+1;
    tr[pos].p00=tr[tr[pos].tr[0]].p00*tr[pos].v0*tr[tr[pos].tr[1]].p00;
    tr[pos].p01=tr[tr[pos].tr[0]].p01*tr[pos].v1*tr[tr[pos].tr[1]].p01;
    tr[pos].p10=tr[tr[pos].tr[1]].p10*tr[pos].v0*tr[tr[pos].tr[0]].p10;
    tr[pos].p11=tr[tr[pos].tr[1]].p11*tr[pos].v1*tr[tr[pos].tr[0]].p11;
  }
  void split(int pos,int k,int&a,int&b){
    if(!pos){
      a=b=0;
      return;
    }
    pushdown(pos);
    if(tr[tr[pos].tr[0]].sz<k)a=pos,split(tr[pos].tr[1],k-tr[tr[pos].tr[0]].sz-1,tr[a].tr[1],b);
    else b=pos,split(tr[pos].tr[0],k,a,tr[b].tr[0]);
    pushup(pos);
  }
  int merge(int a,int b){
    if(!a||!b)return a^b;
    if(tr[a].r<tr[b].r)return pushdown(a),tr[a].tr[1]=merge(tr[a].tr[1],b),pushup(a),a;
    else return pushdown(b),tr[b].tr[0]=merge(a,tr[b].tr[0]),pushup(b),b;
  }
  void append(int&rt,char opt){
    tr[++cnt].sz=1,tr[cnt].r=rand(),tr[cnt].v0=tr[cnt].p00=tr[cnt].p10=opt=='W'?W:E,tr[cnt].v1=tr[cnt].p01=tr[cnt].p11=opt=='W'?E:W,rt=merge(rt,cnt);
  }
  void reverse(int&rt,int l,int r,int a=0,int b=0,int c=0){
    split(rt,l-1,a,b),split(b,r-l+1,b,c),pushdown1(b),rt=merge(merge(a,b),c);
  }
  void flip(int&rt,int l,int r,int a=0,int b=0,int c=0){
    split(rt,l-1,a,b),split(b,r-l+1,b,c),pushdown2(b),rt=merge(merge(a,b),c);
  }
};
FHQtreap<200005>t;
int main(){
  ios::sync_with_stdio(0),cin.tie(0),srand(time(0)),cin>>n>>m>>s;
  for(int i=0;i<n;i++)t.append(rt,s[i]);
  cout<<t.tr[rt].p00.z<<' '<<(t.tr[rt].p00.x+t.tr[rt].p00.z)%mod<<'\n';
  for(int i=1,l,r;i<=m;i++){
    cin>>opt;
    if(opt[0]=='A')cin>>s,t.append(rt,s[0]);
    else if(opt[0]=='F')cin>>l>>r,t.flip(rt,l,r);
    else cin>>l>>r,t.reverse(rt,l,r);
    cout<<t.tr[rt].p00.z<<' '<<(t.tr[rt].p00.x+t.tr[rt].p00.z)%mod<<'\n';
  }
  return 0;
}

Stern-Brocot 树

定义 Stern-Brocot 树的开始是两个分数 \frac 01,\frac 10，可以看作第 0 层。之后的每一层都是在上一层的基础上，每相邻两个分数 \frac ab,\frac cd 中间插入 \frac{a+c}{b+d}。第 k 层称为 k 阶 Stern–Brocot 序列。

假如取出每一层新增的数，会形成一个完全二叉树结构。

Stern-Brocot 树的每一行都是单调递增的，且每个分数都是最简分数（\frac 01,\frac 10 视作 0,\infty）。证明可以归纳，对于两个分数 \frac ab\leq\frac cd，其后代为 \frac{a+c}{b+d}，有 \frac ab\leq\frac{a+c}{b+d}\leq\frac cd。另外，可以归纳出相邻两个分数有 bc-ad=1。对于 \frac ab,\frac{a+c}{b+d}，因为 \frac ab 最简，这个方程的解满足 \frac{a+c}{b+d} 最简。

因此，Stern-Brocot 树可以看成二叉查找树，可以在上面二分。在树上走的过程可以看作有一个三元组 (\frac ab,\frac{a+c}{b+d},\frac cd)，走到 (\frac ab,\frac{2a+c}{2b+d},\frac{a+c}{b+d}) 或 (\frac{a+c}{b+d},\frac{a+2c}{b+2d},\frac cd)，(\frac ab,\frac cd) 也是子树的范围。

值得注意的是，查找一个数，拐弯是次数是 O(\log n) 的，因为连续拐弯会形成斐波那契数列，因此可以将连续的向左或右压成一段。则一个分数可以表示为操作序列 (k_i,opt_i)，其中 opt_i\in\{\texttt L,\texttt R\}，表示 a\gets a+k_ic,b\gets b+k_id 或 c\gets c+k_ia,d\gets d+k_ib。若我们要逼近一个整数，只需每次解不等式求出前进的步数。

假如查找一个分数，有更简单的办法。不妨强制 opt 形如 \texttt{RLRLR}\cdots，在第一步 k_i 可以为 0，其余 k 都为正。假设最后子树范围的端点移动到 \frac pq，设第 i 次移动子树范围的端点的位置是 \frac{p_i}{q_i}，则 开始移动时两端是 \frac{p_{i-1}}{q_{i-1}},\frac{p_{i-2}}{q_{i-2}}，有 \frac{p_i}{q_i}=\frac{k_ip_{i-1}+p_{i-2}}{k_iq_{i-1}+q_{i-2}}。发现这和渐进分数的递推表示相同，设 \frac pq=[a_0,a_1,\cdots,a_n]，则 a 恰好也是在 Stern-Brocot 树上走到 \frac pq 的子节点的路径。若要走到 \frac pq，只需要将 a_n 减一即可。

若我们要求的数是未知的满足某个限制的数，直接二分步长，是 O(\log^2 n) 的。考虑先倍增，当 2^i 不满足限制时，得到上下界 [2^{i-1},2^i]，然后在里面二分或倍增，这样可以保证每次询问次数是 \log 步长。可以证明这样是 O(\log n) 的：考虑相邻两次步长为 u,v，则走完后变为 (\frac{a+vc+uva}{b+vd+uvd},\frac{c+ua}{d+ub})。发现 a,b 增大了 uv 倍。

vector<pair<long long,char> >encode_path(long long p,long long q){
  bool f=1;
  vector<pair<long long,char> >ans;
  while(q){
    if(p>=q)ans.push_back(make_pair(p/q,f?'R':'L'));
    p%=q,swap(p,q),f^=1;
  }
  if(!--ans.back().first)ans.pop_back();
  return ans;
}

Stern-Brocot 树的另一个应用是双曲线下数点。一个经典问题是求 \sum_{i=1}^nd(i)=\sum_{xy\leq n}1，也就是第一象限 y=\frac nx 下方的整点数。这可以做到 O(\sqrt[3]n\log n)。

首先可以用狄利克雷双曲线法转化为 2\sum_{i=1}^{\lfloor\sqrt n\rfloor}\lfloor\frac ni\rfloor-\lfloor\sqrt n\rfloor^2。

考虑用一个凸包拟合这个曲线，令 B=\lfloor\sqrt n\rfloor，初始时有一个点 (\lfloor\frac nB\rfloor+1,B)，每次沿一个向量走，要求走完后是在曲线外面的整点且没有点在向量和曲线中间，将答案加上向量的贡献。

如图，向量 \vec{AB} 的贡献是梯形区域 ABCD（不算 A,B 和线段 BD）的整点数。设从 (x,y) 走到 (x+\Delta x,y-\Delta y)， 贡献是 x\Delta y+\frac{(\Delta x-1)(\Delta +1)}2-\Delta x。

维护一个斜率单调递减的栈，初始为 (1,0),(1,1)。每次先不断沿着栈顶的向量走。当即将走进曲线内，就不断出栈，找到两个向量 \mathbf l,\mathbf r，使得沿 \mathbf l 走能走进曲线，而沿 \mathbf r 走不能。然后在 Stern-Brocot 树上二分出下一步要走的向量。

设 \mathbf{mid} 为 \mathbf l,\mathbf r 在 Stern-Brocot 树上的后代，也就是 \mathbf l+\mathbf r。如果沿 \mathbf{mid} 走不进区间，就令 \mathbf r\gets\mathbf{mid}，同时入栈 \mathbf r；否则 \mathbf l\gets\mathbf{mid}，如果 \mathbf r 的斜率不大于曲线在 x+\mathbf{mid}_x 处的导数，则结束二分。这是因为接下来二分，假如都是向外走，也就是不断令 \mathbf{mid}\gets\mathbf{mid}+\mathbf r。因为 \mathbf r 的斜率已经不超过曲线的斜率，所以不能走到曲线外，\mathbf r 就是下一步走的向量。否则 \mathbf l\gets\mathbf{mid} 接着二分。

当 y\leq n^{\frac 13}，剩下的暴力。可以证明复杂度是 O(\sqrt[3]n\log n)。

推广到 \sum_{xy\leq n}x^py^q 也是可做的。仍然用双曲线法化为 \sum_{i=1}^Bi^pS_q(\lfloor\frac ni\rfloor)+\sum_{i=1}^Bi^qS_p(\lfloor\frac ni\rfloor)，以前一部分为例，从 (x,y) 移动到 (x+\Delta x,y-\Delta y) 的贡献是 -x^py^q+\sum_{i=0}^{\Delta y-1}(y-i)^pS_q(x+\lfloor\frac{i\Delta x}{\Delta y}\rfloor)。自然数幂和是 q+1 次多项式，展开后会变成 O(pq) 个形如 f_{p,q}(x,y)=\sum_{i=0}^{y-1}i^p\lfloor\frac{xi}y\rfloor^q 的问题。

用类欧或万欧是可行的。更好的方式是，整个过程中我们走的向量都由先前的两个向量 \mathbf l,\mathbf r 合并出来，而因为 \mathbf l,\mathbf r 是 Stern-Brocot 树上的邻居，所以 \mathbf l,\mathbf r,\mathbf{mid} 围成的三角形中没有整点。因此可以合并，将 \mathbf{mid} 换成 \mathbf l,\mathbf r 拼起来，有 f_{p,q}(\mathbf{mid})=f_{p,q}(\mathbf l)+\sum_{i=0}^{\mathbf r_y-1}(\mathbf l_y+i)^p(\mathbf l_x+\lfloor\frac{i\mathbf r_x}{\mathbf r_y}\rfloor)^q。同样可以展开，用 \mathbf l,\mathbf r 的信息合并。复杂度 O(\sqrt[3]n\log np^2q^2)。

我们也可以拟合其他的曲线，上述算法的条件是曲线下凸且斜率总不小于 -\frac12。

P5179

考虑在 Stern-Brocot 树上二分，小于下界就向右，大于上界就向左。这个最坏是 O(n) 的。

计算一段连续向右的次数，解不等式 \frac{x+ct}{y+dt}\geq \frac AB，即 t\geq\frac{Ay-Bx}{Bc-Ad}。同理可得向左的次数为 \frac{Cy-Dx}{Da-Cb}。复杂度 O(\log n)。

#include<bits/stdc++.h>
using namespace std;
int t,A,B,C,D;
void solve(int a,int b,int c,int d){
  int x=a+c,y=b+d,temp;
  if(1ll*A*y<1ll*B*x&&1ll*x*D<1ll*y*C){
    cout<<x<<'/'<<y<<'\n';
    return;
  }
  if(1ll*A*y>=1ll*B*x)temp=(1ll*A*y-1ll*B*x)/(1ll*B*c-1ll*A*d),solve(x+temp*c,y+temp*d,c,d);
  else temp=(1ll*C*y-1ll*D*x)/(1ll*D*a-1ll*C*b),solve(a,b,x+temp*a,y+temp*b);
}
int main(){
  while(cin>>A>>B>>C>>D)solve(0,1,1,0);
  return 0;
}

P1797

上文我们 O(\log n) 实现了 ENCODE_PATH 和 DECODE_PATH。LCA 可以取两个分数的公共路径后 DECODE_PATH，ANCESTOR 截取路径前 k 步，RANGE 就是 DECODE_PATH 最后的三元组中的左右两个。

#include<bits/stdc++.h>
using namespace std;
int t,p,q,k,u,v;
char c;
string opt;
vector<pair<int,char> >encode_path(int p,int q){
  bool f=1;
  vector<pair<int,char> >ans;
  while(q){
    if(p>=q)ans.push_back(make_pair(p/q,f?'R':'L'));
    p%=q,swap(p,q),f^=1;
  }
  if(!--ans.back().first)ans.pop_back();
  return ans;
}
pair<pair<int,int>,pair<int,int> >range(vector<pair<int,char> >e){
  int a=0,b=1,c=1,d=0;
  for(int i=0;i<e.size();i++){
    if(e[i].second=='L')c+=a*e[i].first,d+=b*e[i].first;
    else a+=c*e[i].first,b+=d*e[i].first;
  }
  return make_pair(make_pair(a,b),make_pair(c,d));
}
pair<int,int>decode_path(vector<pair<int,char> >e){
  pair<pair<int,int>,pair<int,int> >temp=range(e);
  return make_pair(temp.first.first+temp.second.first,temp.first.second+temp.second.second);
}
pair<int,int>lca(int a,int b,int c,int d){
  vector<pair<int,char> >e,e1=encode_path(a,b),e2=encode_path(c,d);
  for(int i=0;i<min(e1.size(),e2.size());i++){
    if(e1[i]==e2[i])e.push_back(e1[i]);
    else if(e1[i].second!=e2[i].second)break;
    else{
      e.push_back(make_pair(min(e1[i].first,e2[i].first),e1[i].second));
      break;
    }
  }
  return decode_path(e);
}
pair<int,int>ancestor(int p,int q,int k){
  int dep=0;
  vector<pair<int,char> >e=encode_path(p,q);
  for(int i=0;i<e.size();i++)dep+=e[i].first;
  if(k>dep)return make_pair(-1,-1);
  k=dep-k;
  while(k){
    if(k>=e.back().first)k-=e.back().first,e.pop_back();
    else e.back().first-=k,k=0;
  }
  return decode_path(e);
}
int main(){
  cin>>t;
  while(t--){
    cin>>opt;
    if(opt=="ENCODE_PATH"){
      cin>>p>>q;
      vector<pair<int,char> >e=encode_path(p,q);
      cout<<e.size()<<' ';
      for(int i=0;i<e.size();i++)cout<<e[i].second<<' '<<e[i].first<<' ';
      putchar('\n');
    }
    else if(opt=="DECODE_PATH"){
      vector<pair<int,char> >e;
      cin>>k;
      for(int i=0;i<k;i++)cin>>c>>p,e.push_back(make_pair(p,c));
      pair<int,int>temp=decode_path(e);
      cout<<temp.first<<' '<<temp.second<<'\n';
    }
    else if(opt=="LCA"){
      cin>>p>>q>>u>>v;
      pair<int,int>temp=lca(p,q,u,v);
      cout<<temp.first<<' '<<temp.second<<'\n';
    }
    else if(opt=="ANCESTOR"){
      cin>>k>>p>>q;
      pair<int,int>temp=ancestor(p,q,k);
      if(temp.first==-1)puts("-1");
      else cout<<temp.first<<' '<<temp.second<<'\n';
    }
    else if(opt=="RANGE"){
      cin>>p>>q;
      vector<pair<int,char> >e=encode_path(p,q);
      pair<pair<int,int>,pair<int,int> >temp=range(e);
      cout<<temp.first.first<<' '<<temp.first.second<<' '<<temp.second.first<<' '<<temp.second.second<<'\n';
    }
  }
  return 0;
}

SP26073

双曲线下数点模板题。

#include<bits/stdc++.h>
using namespace std;
template<typename T>void out(T x){
  if(x>9)out(x/10);
  putchar(x%10+'0');
}
int t,top;
long long n;
struct vec{
  long long x,y;
}st[1000005];
__int128 solve(long long n,__int128 ans=0){
  long long bn=sqrt(n),bn1=cbrt(n),x=n/bn+1,y=bn;
  top=0,st[++top]=(vec){1,0},st[++top]=(vec){1,1};
  while(1){
    vec l=st[top],r=st[--top];
    while((__int128)(x+l.x)*(y-l.y)>n)ans+=(__int128)x*l.y+(__int128)(l.x-1)*(l.y+1)/2-l.x,x+=l.x,y-=l.y;
    if(y<=bn1)break;
    while((__int128)(x+r.x)*(y-r.y)<=n)l=r,r=st[--top];
    while(1){
      vec mid=(vec){l.x+r.x,l.y+r.y};
      if((__int128)(x+mid.x)*(y-mid.y)>n)st[++top]=mid,r=mid;
      else{
        if((__int128)n*r.x<=(__int128)(x+mid.x)*(x+mid.x)*r.y)break;
        l=mid;
      }
    }
  }
  for(int i=1;i<=y;i++)ans+=n/i;
  return 2*ans-bn*bn;
}
int main(){
  cin>>t;
  while(t--)cin>>n,out(solve(n)),putchar('\n');
  return 0;
}

SP33039

双曲线法转化为 \sum_{i=1}^Bi\lfloor\frac ni\rfloor+\sum_{i=1}^B\frac{\lfloor\frac ni\rfloor(\lfloor\frac ni\rfloor+1)}2-\frac{B^2(B+1)}2。也就是 p=0,q=1 和 p=1,q=0 的问题。

展开走一步的贡献：

\begin{aligned}-y+\sum_{i=0}^{\Delta y-1}(y-i)(x+\lfloor\frac{i\Delta x}{\Delta y}\rfloor)&=-y+xy\Delta y+yf_{0,1}(\mathbf l)-\frac{y\Delta y(\Delta y-1)}2-f_{1,1}(\mathbf l)\\-x+\sum_{i=0}^{\Delta y-1}\frac12(x+\lfloor\frac{i\Delta x}{\Delta y}\rfloor)(x+\lfloor\frac{i\Delta x}{\Delta y}\rfloor+1)&=-x+\frac12(\Delta yx^2+(2x+1)f_{0,1}(\mathbf l)+f_{0,2}(\mathbf l))\end{aligned}

展开 f_{0,1}(\mathbf{mid}),f_{1,1}(\mathbf{mid}),f_{0,2}(\mathbf{mid})：

\begin{aligned}f_{0,1}(\mathbf{mid})&=f_{0,1}(\mathbf l)+\mathbf l_x\mathbf r_y+f_{0,1}(\mathbf r)\\f_{1,1}(\mathbf{mid})&=f_{1,1}(\mathbf l)+\sum_{i=0}^{\mathbf r_y-1}(\mathbf l_y+i)(\mathbf l_x+\lfloor\frac{i\mathbf r_x}{\mathbf r_y}\rfloor)\\&=f_{1,1}(\mathbf l)+\mathbf l_x\mathbf l_y\mathbf r_y+\mathbf l_yf_{0,1}(\mathbf r)+\frac{\mathbf l_x\mathbf r_y(\mathbf r_y-1)}2+f_{1,1}(\mathbf r)\\f_{0,2}(\mathbf{mid})&=f_{0,2}(\mathbf l)+\sum_{i=0}^{\mathbf r_y-1}(\mathbf l_x+\lfloor\frac{i\mathbf r_x}{\mathbf r_y}\rfloor)^2\\&=f_{0,2}(\mathbf l)+\mathbf r_y\mathbf l_x^2+2\mathbf l_xf_{0,1}(\mathbf r)+f_{0,2}(\mathbf r)\end{aligned} 
#include<bits/stdc++.h>
using namespace std;
template<typename T>void out(T x){
  if(x>9)out(x/10);
  putchar(x%10+'0');
}
int t,top;
long long n;
struct vec{
  long long x,y;
  __int128 f01,f11,f02;
}st[1000005];
__int128 solve(long long n){
  long long bn=sqrt(n),bn1=cbrt(n),x=n/bn+1,y=bn;
  __int128 ans=0;
  top=0,st[++top]=(vec){1,0,0,0,0},st[++top]=(vec){1,1,0,0,0};
  while(1){
    vec l=st[top],r=st[--top];
    while((__int128)(x+l.x)*(y-l.y)>n){
      ans+=(__int128)l.y*x*y+y*l.f01-(__int128)l.y*(l.y-1)/2*x-l.f11;
      ans+=((__int128)l.y*x*(x+1)+(2*x+1)*l.f01+l.f02)/2;
      ans-=x+y;
      x+=l.x,y-=l.y;
    }
    if(y<=bn1)break;
    while((__int128)(x+r.x)*(y-r.y)<=n)l=r,r=st[--top];
    while(1){
      vec mid;
      mid.x=l.x+r.x,mid.y=l.y+r.y;
      mid.f01=l.f01+(__int128)l.x*r.y+r.f01;
      mid.f11=l.f11+(__int128)l.x*r.y*(r.y-1)/2+r.f11+(__int128)r.y*l.y*l.x+l.y*r.f01;
      mid.f02=l.f02+(__int128)r.y*l.x*l.x+(__int128)2*l.x*r.f01+r.f02;
      if((__int128)(x+mid.x)*(y-mid.y)>n)r=mid,st[++top]=r;
      else{
        if((__int128)n*r.x<=(__int128)(x+mid.x)*(x+mid.x)*r.y)break;
        l=mid;
      }
    }
  }
  for(int i=1;i<=y;i++)ans+=(__int128)i*(n/i)+(__int128)(n/i)*(n/i+1)/2;
  return ans-(__int128)bn*bn*(bn+1)/2;
}
int main(){
  cin>>t;
  while(t--)cin>>n,out(solve(n)-(__int128)n*(n+1)/2),putchar('\n');
  return 0;
}

P8058

先考虑怎么求一个分数的排名：

\begin{aligned}\sum_{i=1}^n\sum_{j=1}^i[\gcd(i,j)=1][\frac ij\leq\frac pq]&=\sum_{i=1}^n\sum_{j=1}^i[\frac ij\leq\frac pq]\sum_{d\mid i,d\mid j}\mu(d)\\&=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^i[\frac ji\leq\frac pq]\\&=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\lfloor\frac nd\rfloor}\lfloor\frac{pi}q\rfloor\end{aligned}

整除分块，\mu(d) 的块筛用杜教筛算，\sum_{i=1}^{\lfloor\frac nd\rfloor}\lfloor\frac{pi}q\rfloor 用类欧算，O(\sqrt n\log n)。有个小优化，每次算之前先预处理一段 \sum_{i=1}^{\lfloor\frac nd\rfloor}\lfloor\frac{pi}q\rfloor，变为 O(\sqrt{n\log n})。

外层我们需要二分这个分数。考虑 Stern-Brocot 树上倍增，向左走时倍增第一个 \frac{a+c}{b+d} 的排名 \leq k 的位置，向右走时倍增第一个 \frac{a+c}{b+d} 的排名 \geq k 的位置。总复杂度 O(\sqrt n\log^{1.5}n)。

#include<bits/stdc++.h>
using namespace std;
int bn,bn1,prime[1000005],mu[1000005];
long long n,k,pre[173210];
bool vis[1000005];
unordered_map<long long,int>sum;
int mu_init(int n,int prime[],int mu[],bool vis[],int cnt=0){
  mu[1]=1;
  for(int i=2;i<=n;i++)vis[i]=1;
  for(int i=2;i<=n;i++){
    if(vis[i])prime[++cnt]=i,mu[i]=-1;
    for(int j=1,v,p;i*prime[j]<=n;j++){
      p=prime[j],v=i*p,vis[v]=0;
      if(i%p==0){
        mu[v]=0;
        break;
      }
      else mu[v]=-mu[i];
    }
  }
  return cnt;
}
long long getsum(long long x){
  if(x<=bn)return mu[x];
  if(sum.count(x))return sum[x];
  long long ans=1;
  for(long long l=2,r;l<=x;l=r+1)r=x/(x/l),ans-=(r-l+1)*getsum(x/l);
  return sum[x]=ans;
}
long long solve(long long n,long long a,long long b,long long c){
  if(!a)return (n+1)*(b/c);
  else if(a>=c||b>=c)return solve(n,a%c,b%c,c)+n*(n+1)/2*(a/c)+(n+1)*(b/c);
  else return n*((a*n+b)/c)-solve((a*n+b)/c-1,c,c-b-1,a);
}
long long calc(long long p,long long q,long long ans=0){
  for(int i=1;i<=bn1;i++)pre[i]=pre[i-1]+p*i/q;
  for(long long l=1,r;l<=n;l=r+1)r=n/(n/l),ans+=(getsum(r)-getsum(l-1))*(n/l<=bn1?pre[n/l]:solve(n/l,p,0,q));
  return ans;
}
int main(){
  cin>>n>>k,bn=pow(n,2.0/3),bn1=sqrt(n*log2(n)),mu_init(bn,prime,mu,vis);
  for(int i=1;i<=bn;i++)mu[i]+=mu[i-1];
  bool f=1;
  long long a=0,b=1,c=1,d=1;
  while(1){
    if(f){
      int p=1;
      while(calc(a+c*p+c,b+d*p+d)<k)p<<=1;
      p>>=1,a+=c*p,b+=d*p;
      for(p>>=1;p>=1;p>>=1)if(calc(a+c*p+c,b+d*p+d)<k)a+=c*p,b+=d*p;
      if(calc(a+c,b+d)<k)a+=c,b+=d;
    }
    else{
      int p=1;
      while(calc(a+c+a*p,b+d+b*p)>k)p<<=1;
      p>>=1,c+=a*p,d+=b*p;
      for(p>>=1;p>=1;p>>=1)if(calc(a+c+a*p,b+d+b*p)>k)c+=a*p,d+=b*p;
      if(calc(a+c,b+d)>k)c+=a,d+=b;
    }
    if(calc(a+c,b+d)==k)break;
    f^=1;
  }
  return cout<<a+c<<' '<<b+d<<'\n',0;
}

P1676

我们用 Stern-Brocot 树去拟合正方形减掉圆形的一个角，即 y=n-\sqrt{2nx-x^2}。

我们需要先分成斜率 \geq-\frac12 的两部分。用直线 y=x 与曲线交，最近的整点是 B=\lfloor(1-\frac{\sqrt2}2)n\rfloor，将整点分成 x\leq B,y\leq B，最后加上 B^2。此时这半条曲线可以从 (\lceil n-\sqrt{2nB+B^2}\rceil,B)（横坐标是最小的 x 不严格在曲线内部）开始拟合。过程基本不变，f'(x)=-\frac{n-x}{\sqrt{2nx-x^2}}。最后取一个阈值开始跑暴力。

这个复杂度未知，但可以说明不超过 O(n^\frac23)，也就是平面凸包点数。只考虑右下凸包。放缩一下，忽略约分后相同的情况，贪心取分子分母和最小的几个作为组成凸包的向量，必要条件是终点在 [0,n]^2 之内。贪心地取 (1,1),(1,2),(2,1),(1,3),(2,2),(3,1)\cdots，容易计算前 k 项的和是 O(k^\frac 32) 的。而终点的坐标和不超过 2n，因此 k=O(n^\frac 23)。

#include<bits/stdc++.h>
using namespace std;
template<typename T>void out(T x){
  if(x>9)out(x/10);
  putchar(x%10+'0');
}
int top;
long long n;
struct vec{
  long long x,y;
}st[10000005];
bool checkin(long long x,long long y,long long n){
  return y<0||(__int128)x*(2*n-x)<(__int128)(n-y)*(n-y);
}
long double der(long long x,long long n){
  return (x-n)/sqrtl((__int128)x*(2*n-x));
}
__int128 solve(long long n,__int128 ans=0){
  long long bn=(1-M_SQRT1_2)*n,bn1=sqrt(n),x=ceil(n-sqrtl((__int128)bn*(2*n-bn))),y=bn;
  top=0,st[++top]=(vec){1,0},st[++top]=(vec){1,1};
  while(y){
    vec l=st[top],r=st[--top];
    while(!checkin(x+l.x,y-l.y,n))ans+=(__int128)x*l.y+(__int128)(l.x-1)*(l.y+1)/2-l.x,x+=l.x,y-=l.y;
    if(y<=bn1)break;
    while(checkin(x+r.x,y-r.y,n))l=r,r=st[--top];
    while(1){
      vec mid=(vec){l.x+r.x,l.y+r.y};
      if(!checkin(x+mid.x,y-mid.y,n))st[++top]=mid,r=mid;
      else{
        if(r.y>-r.x*der(x+mid.x,n))break;
        l=mid;
      }
    }
  }
  for(int i=1;i<=y;i++)ans+=n-(__int128)sqrtl((__int128)i*(2*n-i))-1;
  return 2*ans-(__int128)bn*bn;
}
int main(){
  return cin>>n,out((__int128)4*n*n-4*solve(n)-4*n+5),0;
}

Farey 序列

定义 Farey 序列 \mathcal{F}_n 为分母不超过 n 的最简真分数的升序排列。特殊地，认为 \frac 0 1,\frac 1 1 也是最简真分数。比如 \mathcal{F}_5=\left\{\frac 0 1,\frac 1 5,\frac 1 4,\frac 1 3,\frac 2 5,\frac 1 2,\frac 3 5,\frac 2 3,\frac 3 4,\frac 4 5,\frac 1 1\right\}。

考虑 Farey 序列和 Stern-Brocot 树的关系，在每次插入时，只需把产生的分母大于 n 的分数删掉，就得到 Farey 序列，因此也满足单调性和最简性。