连分数

设实数 \(x=[a_0,a_1,a_2,\cdots]\)。那么，成立如下性质：

1. 对于任意 \(k\in\mathbf Z\)，有 \(x+k=[a_0+k,a_1,a_2,\cdots]\)；
2. 对实数 \(x>1\)，有 \(a_0>0\)，且它的倒数 \(x^{-1}=[0,a_0,a_1,a_2,\cdots]\)。

渐近分数 \(x_k\) 的分子和分母可以看作 \(a_0, a_1, \cdots, a_k\) 的多元多项式：

根据渐近分数的定义，有

和上式比较，不妨设 \(Q_k(a_0, \cdots, a_k) = P_{k-1}(a_1, \cdots, a_k)\)，就可以将渐近分数写作

且多项式 \(P_k\) 有递推关系

因为

所以，递推的起点是

如果设

可以验证对于 \(k=0,1\) 也成立上述递推关系。这相当于规定了形式分数 \(r_{-1}=\dfrac{1}{0}\) 和 \(r_{-2}=\dfrac{0}{1}\)。

满足上述递推关系的多项式列 \(P_k\) 称为 连项式²（continuant）。它可以写作行列式的形式：

这是一个 三对角矩阵 的行列式，从左上角开始展开，可以验证它具有上面的递推关系和初值条件。反过来，从右下角开始展开，则又能得到递推关系

这就是所要求证的。

设实数 \(x=[a_0,a_1,a_2,\cdots]\) 的第 \(k\) 个渐近分数是 \(\dfrac{p_k}{q_k}\)，则相邻两个渐近分数的分子和分母的比值分别为

如果 \(a_0=0\)，则第一个连分数应当理解为在倒数第二项处截断，即 \([a_k,a_{k-1},\cdots,a_2]\)。

在 \(p_k\) 和 \(q_k\) 的递推关系中，左右两侧分别同除以 \(p_{k-1}\) 和 \(q_{k-1}\)，就得到

迭代这两个式子，就可以得到两个连分数。再代入初始值 \(\dfrac{p_0}{p_{-1}}=a_0\) 和 \(\dfrac{q_1}{q_0}=a_1\) 即可。至于 \(a_0=0\) 的情形，将得到的连分数理解为形式表达式，则它的余项

因而可以直接略去最后两项。如果需要严格的证明，只要注意到这个式子可以看作 \(a_0\rightarrow 0\) 的极限即可。

实数 \(x>0\) 的渐近分数的倒数是 \(x^{-1}\) 的渐近分数。

不妨设 \(x>1\) 且有连分数表示 \([a_0,a_1,a_2,\cdots]\)，则 \(x^{-1}\) 的连分数表示是 \([0,a_0,a_1,a_2,\cdots]\)。它们的渐近分数可以从递推关系中求得。而且，对于 \(x\)，有初值条件 \(x_{-2}=\dfrac{0}{1}\) 和 \(x_{-1}=\dfrac{1}{0}\)；对于 \(y=x^{-1}\)，有初值条件 \(y_{-1}=\dfrac{1}{0}\) 和 \(y_{0}=\dfrac{0}{1}\)。因此，有 \(x_{-2}=(y_{-1})^{-1}\) 和 \(x_{-1}=(y_0)^{-1}\)。根据递推关系，可以得到 \(x_k=y_{k+1}^{-1}\)。这就说明，\(x\) 的渐近分数的倒数是 \(y=x^{-1}\) 的渐近分数。对于 \(0<x\le 1\) 的情形，也可以做类似讨论。

设 \(x_k=\dfrac{p_k}{q_k}\) 是实数 \(x\) 的第 \(k\) 个渐近分数。那么，有

因此，相邻两项的渐近分数的差分是

根据递推关系，有

这就是 \(p_{k+1}q_k − p_kq_{k+1} = (−1)^k\)。对两边同除以 \(q_{k+1}q_k\)，就得到关于 \(x_{k+1}-x_k\) 的结论。

设 \(x_k=\dfrac{p_k}{q_k}\neq x\) 是实数 \(x\) 的第 \(k\) 个渐近分数。那么，有

其中，\(r_{k+1}\) 是实数 \(x\) 的第 \(k+1\) 个余项。进而，有

因为 \(x=[a_0,a_1,\cdots,a_k,r_{k+1}]\)，而且对形式连分数也成立渐近分数的差分公式，所以有

其中，\(r_{k+1}q_k+q_{k-1}\) 就是按照递推公式得到的这个形式连分数的第 \(k+1\) 个渐近分数的分母。

要完成随后的不等式估计，只需要注意到当 \(x_k\neq x\) 时，总成立

所以有

因此，有不等式

要得到外侧的放缩，再注意到 \(q_k\le q_{k+1}\) 就可以了。

对于任何实数 \(x\)，且渐近分数 \(x_k=\dfrac{p_k}{q_k}\) 的分子和分母由递推公式给出，则 \(\dfrac{p_k}{q_k}\) 是既约分数，即 \(\gcd(p_k,q_k)=1\)。

对差分公式应用 裴蜀定理 即可。

虽然这个问题通常是用 扩展欧几里得算法 解决的，但是同样可以通过连分数求解。

设 \(\dfrac{A}{B}=[a_0, a_1, \cdots, a_k]\)。上面证明了 \(p_k q_{k-1} - p_{k-1} q_k = (-1)^{k-1}\)。将 \(p_k\) 和 \(q_k\) 替换为 \(A\) 和 \(B\)，得到

其中，\(g = \gcd(A, B)\)。如果 \(g\) 整除 \(C\)，则一组解为 \(x = (-1)^{k-1}\dfrac{C}{g} q_{k-1}\) 和 \(y = (-1)^{k}\dfrac{C}{g} p_{k-1}\)；否则无解。

1
2
3
4
5
6
7
8
9

// Return (x,y) such that Ax+By=C.
// Assume that such (x,y) exists.
auto dio(int A, int B, int C) {
  std::vector<int> p, q;
  std::tie(p, q) = convergents(fraction(A, B));
  C /= A / p.back();
  int t = p.size() % 2 ? -1 : 1;
  return std::make_pair(t * C * q.end()[-2], -t * C * p.end()[-2]);
}

1
2
3
4
5
6
7

# Return (x, y) such that Ax+By=C.
# Assume that such (x, y) exists.
def dio(A, B, C):
    p, q = convergents(fraction(A, B))
    C //= A // p[-1]  # divide by gcd(A, B)
    t = (-1) if len(p) % 2 else 1
    return t * C * q[-2], -t * C * p[-2]

对于无理数 \(x\)，存在无穷多个既约分数 \(\dfrac{p}{q}\) 使得

成立。

根据渐近分数的误差估计，对于无理数 \(x\) 的第 \(k\) 个渐近分数 \(x_k=\dfrac{p_k}{q_k}\)，有

检查误差公式的证明即可知，对于任一无理数 \(x\)，取等条件并不成立。因此，它的所有渐近分数的分子和分母都满足要求。

对于无理数 \(x\)，存在无穷多个既约分数 \(\dfrac{p}{q}\) 使得

成立，而且不等式右侧的 \(\sqrt{5}\) 不能换成更大的实数。

Borel 实际上证明了，无理数 \(x\) 的连续三个渐近分数中，必然有至少一个满足上述条件。因为渐近分数无穷多，且都是既约分数，那么 Hurwitz 定理的第一部分就必然成立。

反证法。不妨设存在无理数 \(x\) 和它的渐近分数 \(x_{k-1},x_k,x_{k+1}\)，使得

成立。因为相邻的渐近分数必然位于 \(x\) 两侧，所以由差分公式知

它可以写成关于商 \(\dfrac{q_{k}}{q_{k-1}}\) 的不等式

因为左侧是有理数，右侧是无理数，等号必然无法取得。又因为 \(q_k\ge q_{k-1}\)，所以可以解得

同理，可以证明

但是根据递推公式，并结合两式可知，

这与简单连分数的定义矛盾。所以，Borel 的结论成立。

要说明这样得到的界是最好的，只需要找到 \(x\) 使得对于任何 \(C>\sqrt{5}\)，都只存在有限多个既约分数 \(\dfrac{p}{q}\) 使得不等式

成立。下面证明 \(\phi=\dfrac{\sqrt{5}+1}{2}\) 就是这样的 \(x\)。³

设 \(\phi'=\dfrac{-\sqrt{5}+1}{2}\) 是 \(\phi\) 的共轭根。它们都是方程 \(x^2-x-1=0\) 的根。因而，对任意实数 \(x\)，都有

代入既约分数 \(\dfrac{p}{q}\)，就得到

对于 \(C>\sqrt{5}\)，可以直接解出 \(q<\sqrt{C(C-\sqrt{5})}\)，因而不可能存在无穷多组解满足上述不等式。

所有的第一类最佳逼近都是中间分数。

因为 \(a_0\le x\le a_0+1\)，所以第一类最佳逼近必然位于 \(x_{1,0}=a_0\) 与 \(x_{0,1}=a_0+1\) 之间。所有中间分数从小到大可以排列成

同阶的中间分数是连续出现的，而不同阶的中间分数之间又没有间隔。这意味着，任何位于 \(x_{1,0}=a_0\) 与 \(x_{0,1}=a_0+1\) 之间的有理数 \(\dfrac{p}{q}\) 必然落在两个同阶的中间分数 \(x_{k,t}\) 和 \(x_{k,t+1}\) 之间。不妨设它不是中间分数且小于 \(x\)，因而有

此时，一方面有

另一方面有

因此，必然有

也就是说，有理数 \(\dfrac{p}{q}\) 的分母一定大于 \(x_{k,t+1}\) 的分母，但是它并不是更好的逼近：

因此，它不可能是第一类最佳逼近。这就说明，不是中间分数，就不是第一类最佳逼近；亦即所有第一类最佳逼近都是中间分数。

所有渐近分数都是第一类最佳逼近。除此之外，设 \(0<t<a_{k+1}\)，则中间分数 \(x_{k,t}\) 是第一类最佳逼近，当且仅当 \(t>\dfrac{a_{k+1}}{2}\) 或者 \(t=\dfrac{a_{k+1}}{2}\) 且 \(r_{k+2}>\dfrac{q_k}{q_{k-1}}\)。

下面会证明，渐近分数都是第二类最佳逼近，因而必然是第一类最佳逼近。关键在于那些不是渐近分数的中间分数。

前文已经说过，中间分数 \(x_{k,t}\) 的分母位于 \(x_k\) 和 \(x_{k+1}\) 之间，且随着 \(t\) 的增加逐渐增大，但是 \(x_{k,t}\) 却逐渐接近 \(x_{k+1}\)，因而更接近 \(x\)。以 \(x_{k,t}<x\) 为例，它与相邻的中间分数的相对位置关系满足：

因为 \(x_k\) 的分母小于 \(x_{k,t}\)，所以 \(x_{k,t}\) 成为第一类最佳逼近的必要条件就是，它比 \(x_k\) 更接近 \(x\)。这也是充分条件，因为作为渐近分数，没有比 \(x_k\) 分母更小但是距离更近的了，而那些比 \(x_k\) 分母还要大的中间分数，必然与 \(x_{k,t}\) 同阶，但是分母更小，就必然距离 \(x\) 更远。对于渐近分数和中间分数的误差，经计算可知

其中用到了 \(r_{k+1}\ge a_{k+1}>t\)。因此，\(x_{k,t}\) 比 \(x_k\) 更接近 \(x\)，成为第一类最佳逼近，当且仅当

此时，有三种可能的情况：

1. 如果 \(t<\dfrac{a_{t+1}}{2}\)，那么 \(2t<a_{t+1}\)。因为两侧都是整数，所以 \(2t\le a_{k+1}-1\)，故而 \(2t+\dfrac{q_{k-1}}{q_k}\le 2t+1\le a_{k+1}\le r_{t+1}\)。此时，\(x_{k,t}\) 不是第一类最佳逼近；
2. 如果 \(t>\dfrac{a_{t+1}}{2}\)，那么 \(2t>a_{t+1}\)。因为两侧都是整数，所以 \(2t\ge a_{t+1}+1>r_{t+1}\)。此时，\(x_{k,t}\) 是第一类最佳逼近；
3. 如果 \(a_{t+1}\) 是偶数，还有第三种情况，即 \(t=\dfrac{a_{t+1}}{2}\)。上述条件等价于 \(\dfrac{1}{r_{k+1}}=r_{k+1}-a_{k+1}<\dfrac{q_{k-1}}{q_k}\)，亦即 \(r_{k+2}>\dfrac{q_k}{q_{k-1}}\)。

所有的第二类最佳逼近一定是渐近分数，所有的渐近分数也一定是第二类最佳逼近。

要证明第一部分，因为第二类最佳逼近也一定是第一类最佳逼近，所以只需要证明不是渐近分数的中间分数不能成为第二类最佳逼近就可以了。为此，设 \(x_{k,t}=\dfrac{p}{q}\) 是中间分数但是不是渐近分数，那么，设 \(x_{k,t}<x\)，有

因为 \(x_{k,t}\) 与 \(x\) 之间的误差

并利用渐近分数的误差估计，所以总是有

即 \(x_{k,t}\) 的逼近程度不优于分母更小的 \(x_k\) 逼近程度，所以不可能是第二类最佳逼近。

反过来，要证明第二部分，即每个渐近分数 \(x_k=\dfrac{p_k}{q_k}\) 都是第二类最佳逼近。这就是要说明，对于所有分数 \(\dfrac{p}{q}\) 且 \(q\le q_k\)，都有 \(|q_kx-p_k|<|qx-p|\)。不考虑半奇数的情形，则可以假定 \(k>0\)。首先，根据渐近分数逼近实数的误差估计，有

不等式全部成立等号，当且仅当 \(a_{k+1}=1\) 且是连分数的末项。不考虑这样的情形，那么 \(x_{k-1}=\dfrac{p_{k-1}}{q_{k-1}}\) 严格劣于 \(x_k=\dfrac{p_k}{q_k}\)。

任取一分数 \(\dfrac{p}{q}\neq x_k\) 且 \(0<q\le q_k\)，因为有差分公式 \(p_{k}q_{k-1} − p_{k-1}q_{k} = (−1)^{k-1}\)，所以由 Cramer 法则可知，线性方程组

必然存在唯一的整数解 \((\lambda,\mu)\)。如果 \(\lambda\mu>0\)，那么 \(q>|\lambda|q_k\ge q_k\)，矛盾。否则，\(\lambda\mu\le 0\)，即 \(\lambda\) 和 \(\mu\) 异号，那么因为 \(q_{k-1}x-p_{k-1}\) 和 \(q_kx-p_k\) 也异号，就有 \(\lambda(q_{k-1}x-p_{k-1})\) 和 \(\mu(q_kx-p_k)\) 同号，故而

最后的不等号是严格的，因为 \(x_{k-1}\) 严格劣于 \(x_k\)，且 \(\dfrac{p}{q}\neq x_k\)。这就说明，\(x_k\) 是第二类最佳逼近。

对于实数 \(x\) 与分数 \(\dfrac{p}{q}\)，如果有

那么 \(\dfrac{p}{q}\) 一定是 \(x\) 的渐近分数。

设 \(\epsilon\in\{-1,1\}\) 和 \(\theta\in(0,1/2)\) 是使得

成立的常数。将有理数 \(\dfrac{p}{q}\) 展开成连分数 \([a_0,a_1,\cdots,a_n]\)。此处，有理数有两种连分数表示，其中的 \(n\) 恰好相差一，所以可以取连分数表示使得 \((-1)^n=\epsilon\)，并记这个连分数表示的渐近分数为 \(\dfrac{p_k}{q_k}\)。设实数 \(\omega\) 满足

那么，必然有

故而，有

这说明

将 \(\omega\) 也展成连分数 \([b_0,b_1,\cdots]\)，则

这是合法的简单连分数，所以 \(\dfrac{p}{q}\) 就是 \(x\) 的渐近分数。

这个证明实际说明 \(\dfrac{p}{q}\) 成为渐近分数的充分必要条件是上述证明中的 \(\omega>1\)，这正是 Legendre 判别法的原始形式。

实数 \(x\) 的相邻两个渐近分数中至少有一个满足

假设不然。存在实数 \(x\) 有两个相邻的渐近分数 \(x_{k-1}\) 和 \(x_k\) 满足

因为 \(x\) 位于 \(x_{k-1}\) 和 \(x_k\) 之间，所以

这说明 \(q_k=q_{k+1}\)。故而必然有 \(k=0\) 且 \(a_1=1\)。此时前两个渐近分数是 \(x_0=a_0\) 和 \(x_1=a_0+1\)。故而命题的唯一反例是半奇数，按照前文的说明，本文不考虑这种情形。

• 每个分数 \(\nu=\dfrac{p}{q}\) 都对应着第一象限内的一个整点 \(\vec\nu=(q,p)\)，分数的大小对应着与原点连线的斜率。

• 直线 \(y=\xi x\) 的方向向量是 \(\vec\xi=(1,\xi)\)。利用 叉积 \((x_1,y_1)\times(x_2,y_2)=x_1y_2-x_2y_1\) 的概念，可以通过 \(\vec\xi\times\vec\nu=p-q\xi\) 的正负判断点在直线上方还是下方。因而，在直线上方的点就对应着大于等于 \(\xi\) 的分数，在直线下方的点就对应着小于等于 \(\xi\) 的分数。叉积的绝对值 \(|\vec\xi\times\vec\nu|\) 正比于点 \(\vec\nu\) 与直线 \(y=\xi x\) 的距离

  \[ \dfrac{|p-qx|}{\sqrt{1+\xi^2}}, \]

  对应着分数 \(\nu\) 对实数 \(\xi\) 的逼近程度。

• 将渐近分数 \(\xi_k=\dfrac{p_k}{q_k}\) 对应的点记作 \(\vec\xi_k=(p_k,q_k)\)，则递推公式就可以写作

  \[ \vec\xi_k = a_k\vec\xi_{k-1} + \vec\xi_{k-2}. \]

  递归的起点是 \(\xi_{-2} = (1,0)\) 和 \(\xi_{-1} = (0,1)\)。

• 对于整数 \(t\)，如果 \(0\le t\le a_k\)，那么点

  \[ \vec\xi_{k-1,t} = t\vec\xi_{k-1} + \vec\xi_{k-2} \]

  就落在连结点 \(\vec\xi_{k-2}\) 和点 \(\vec\xi_k\) 的线段上。它们对应着中间分数 \(\xi_{k-1,t}\)。

• 利用几何的方法可以构造出所有的渐近分数和中间分数。从点 \(\vec\xi_{-2}=(1,0)\) 和点 \(\vec\xi_{-1}=(0,1)\) 开始，两个点位于直线 \(y=\xi x\) 两侧，这意味着 \(\vec\xi\times\vec\xi_{-2}\) 和 \(\vec\xi\times\vec\xi_{-1}\) 符号相反。将 \(\vec\xi_{-1}\) 按照向量的加法添加到 \(\vec\xi_{-2}\) 上，直到无法继续添加而不穿过直线 \(y=\xi x\) 为止，将结果记作 \(\vec\xi_0\)，此时仍与 \(\vec\xi_{-1}\) 不同侧。再将 \(\vec\xi_0\) 添加到 \(\vec\xi_{-1}\) 上，直到无法继续添加而不穿过直线 \(y=\xi x\) 为止，将结果记作 \(\vec\xi_1\)，此时仍与 \(\vec\xi_0\) 不同侧。这个过程可以一直持续到无穷，除非在有限步内某个 \(\vec\xi_n\) 就恰好落在直线 \(y=\xi x\) 上。后者意味着向量 \(\vec\xi\) 与 \(\vec\xi_n\) 共线，即 \(\xi=\dfrac{p_n}{q_n}\) 为有理点。这个过程就可以得到前面示意图中的图形。Boris Delaunay 将这个过程形象地称为鼻子拉伸算法（nose-streching algorithm）⁶。

• 如果需要快速计算每一步将 \(\vec\xi_{k-1}\) 添加到 \(\vec\xi_{k-2}\) 需要的次数，可以借助叉积。因为 \(\vec\xi\times\vec\xi_{k-1}\) 与 \(\vec\xi\times\vec\xi_{k-2}\) 符号相反，所以如果记 \(\vec\xi_{k-1,t}=t\vec\xi_{k-1}+\vec\xi_{k-2}\) 为向 \(\vec\xi_{k-2}\) 添加 \(t\) 次 \(\vec\xi_{k-1}\) 得到的结果，则 \(\vec\xi\times\vec\xi_{k-1,t}=t(\vec\xi\times\vec\xi_{k-1})+(\vec\xi\times\vec\xi_{k-2})\) 不改变符号，就意味着没有穿过直线。在不变号之前，\(\vec\xi\times\vec\xi_{k-1,t}\) 的绝对值会逐渐下降。记

  \[ r_{k} = \left|\dfrac{\vec\xi\times\vec\xi_{k-2}}{\vec\xi\times\vec\xi_{k-1}}\right| = -\dfrac{\vec\xi\times\vec\xi_{k-2}}{\vec\xi\times\vec\xi_{k-1}}. \]

  那么，最大可以下降的次数就是

  \[ a_{k} = \lfloor r_{k}\rfloor = \left\lfloor\left|\dfrac{q_{k-1}\xi-p_{k-1}}{q_{k-2}\xi-p_{k-2}}\right|\right\rfloor. \]

  这就是连分数展开的第 \(k\) 项。而且，\(r_k\) 就是连分数展开的余项，它满足关系式：

  \[ r_k = -\dfrac{q_{k-1}\xi-p_{k-1}}{q_{k-2}\xi-p_{k-2}} \iff \xi = \dfrac{p_{k-1}r_k + p_{k-2}}{q_{k-1}r_k+q_{k-2}}. \]

  这就是连分数关系式 \(\xi = [a_0,a_1,\cdots,a_{k-1},r_k]\)。

• 因为每次添加向量造成的 \(\vec\xi\times\vec\xi_{k-1,t}\) 的变化的步长都是 \(|\vec\xi\times\vec\xi_{k-1}|\)，所以最后剩余的距离 \(|\vec\xi\times\vec\xi_k|\) 必然严格小于 \(|\vec\xi\times\vec\xi_{k-1}|\)。这说明，渐近分数的逼近程度（由 \(|qx-p|\) 衡量）是随着 \(k\) 的增加严格更优的。

• 利用叉积的运算法则，有

  \[ \vec\xi_{k}\times\vec\xi_{k+1} = \vec\xi_{k}\times(a_{k+1}\vec\xi_k+\vec\xi_{k-1}) = \vec\xi_{k}\times\vec\xi_{k-1} = -\vec\xi_{k-1}\times\vec\xi_{k}. \]

  归纳可知

  \[ \vec\xi_{k}\times\vec\xi_{k+1} = (-1)^{k+2}\vec\xi_{k-2}\times\vec\xi_{k-1} = (-1)^{k}. \]

  这就是渐近分数的差分公式 \(p_{k+1}q_k-p_kq_{k+1}=(-1)^k\)。

• 上下两个凸包之间的面积可以剖分成若干个（可能是无穷多个）三角形，其中每个三角形的顶点分别是 \(\vec\xi_{k-2}\)、\(\vec\xi_k\) 和 \(\vec 0\)。这样的三角形的面积是

  \[ \dfrac12|\vec\xi_{k-2}\times\vec\xi_k| = \dfrac12|\vec\xi_{k-2}\times(a_k\vec\xi_{k-1}+\vec\xi_{k-2})| = \dfrac{a_k}{2}|\vec\xi_{k-2}\times\vec\xi_{k-1}| = \dfrac{a_k}{2}. \]

  根据 Pick 定理，这意味着如果设 \(I\) 和 \(B\) 分别为三角形内部和边界上的整点个数，则

  \[ I + \dfrac{B}{2} - 1 = \dfrac{a_k}{2}. \]

  又已知三角形边界上已经有了 \(\{\vec 0\}\cup\{\vec\xi_{k-1,t}:0\le t\le a_k\}\) 这共计 \(a_k+2\) 个整点。这说明，就一定有 \(I=0\) 和 \(B=a_k+2\)。因而，三角形的边上没有更多的整点，三角形内部也没有整点。也就是说，\(q_k\) 和 \(p_k\) 是既约的，中间分数是连结 \(\vec\xi_{k-2}\) 和 \(\vec\xi_k\) 的边上的全部整点，且第一象限的所有整点都在上下两个凸包内。

首先将两个连分数表示都转化成末尾是 \(1\) 的形式。不妨设题目所给的已经是这样形式的连分数，即 \(\alpha_n=\beta_m=1\)。因为偶数位置（下标从 \(0\) 开始）是向右移动的步数，奇数位置是向左移动的步数，所以，\(\alpha<\beta\) 当且仅当按照 字典序 比较时，有

相较于连分数表示，交替地添加正负号，删去末尾的 \(1\)，并且长度不足的位置用 \(0\) 补齐。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

// Expand [..., n] to [..., n-1, 1] if needed.
void expand(std::vector<int>& a) {
  if (a.size() == 1 || a.back() > 1) {
    --a.back();
    a.push_back(1);
  }
}

// Check if a is smaller than b.
bool less_than(std::vector<int> a, std::vector<int> b) {
  expand(a);
  expand(b);
  for (int i = 0; i < a.size() - 1 || i < b.size() - 1; ++i) {
    int d = (i < a.size() - 1 ? a[i] : 0) - (i < b.size() - 1 ? b[i] : 0);
    if (i & 1) d = -d;
    if (d < 0) {
      return true;
    } else if (d > 0) {
      return false;
    }
  }
  return false;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

# Expand [..., n] to [..., n-1, 1] if needed.
def expand(a):
    if a[-1] != 1 or len(a) == 1:
        a[-1] -= 1
        a.append(1)
    return a


# Check if a is smaller than b.
def less_than(a, b):
    a = expand(a)
    b = expand(b)
    a = [(-1) ** i * a[i] for i in range(len(a))]
    b = [(-1) ** i * b[i] for i in range(len(b))]
    return a < b

因为 Stern–Brocot 树既是 \([0,\infty]\) 中的分数的二叉搜索树，又是二元组 \((q,p)\) 的 笛卡尔树，所以题意几乎可以转化为求 Stern–Brocot 树上两个点的 LCA（最近公共祖先）。但是，LCA 只能处理闭区间内的情形，LCA 可能是端点本身。为了避免额外的讨论，可以首先构造出 \(\dfrac{p_0}{q_0}+\varepsilon\) 和 \(\dfrac{p_1}{q_1}-\varepsilon\)，再计算 LCA。在已经通过连分数计算出根到节点的路径的情况下，LCA 只要取最长的公共路径即可。

要构造出 \(x\pm\varepsilon\)，只需要在节点 \(x\) 处首先向右（左）移动一次，再向左（右）移动 \(\infty\) 次即可。转化成连分数的语言，对于分数 \(x=[a_0,a_1,\cdots,a_{n-1},1]\)，可以知道 \(x\pm\varepsilon\) 必然是 \([a_0,a_1,\cdots,a_{n-1}+1,\infty]\) 和 \([a_0,a_1,\cdots,a_{n-1},1,\infty]\)，因而只需要比较这两个连分数，将较大（小）的定义为 \(x\pm\varepsilon\)。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

// Get X +- EPSILON.
auto pm_eps(std::vector<int> a) {
  constexpr int inf = 0x3f3f3f3f;
  // Deal with empty continued fraction for 1/0.
  if (a.empty()) {
    a.emplace_back(inf);
  }
  auto b = a;
  expand(b);
  a.emplace_back(inf);
  b.emplace_back(inf);
  return less_than(a, b) ? std::make_pair(a, b) : std::make_pair(b, a);
}

// Find the lexicographically smallest (q, p)
//   such that p0/q0 < p/q < p1/q1.
auto middle(int p0, int q0, int p1, int q1) {
  auto a0 = pm_eps(fraction(p0, q0)).second;
  auto a1 = pm_eps(fraction(p1, q1)).first;
  std::vector<int> a;
  for (int i = 0; i < a0.size() || i < a1.size(); ++i) {
    if (a0[i] == a1[i]) {
      a.emplace_back(a0[i]);
    } else {
      a.emplace_back(std::min(a0[i], a1[i]) + 1);
      break;
    }
  }
  auto pq = convergents(a);
  return std::make_pair(pq.first.back(), pq.second.back());
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

# Get X +- EPSILON.
def pm_eps(a):
    # Deal with empty continued fraction for 1/0.
    if not a:
        a.append(float("inf"))
    b = expand(a.copy())
    a.append(float("inf"))
    b.append(float("inf"))
    return (a, b) if less_than(a, b) else (b, a)


# Find the lexicographically smallest (q, p)
#   such that p0/q0 < p/q < p1/q1.
def middle(p0, q0, p1, q1):
    a0 = pm_eps(fraction(p0, q0))[1]
    a1 = pm_eps(fraction(p1, q1))[0]
    a = []
    for i in range(min(len(a0), len(a1))):
        if a0[i] == a1[i]:
            a.append(a0[i])
        else:
            a.append(int(min(a0[i], a1[i])) + 1)
            break
    p, q = convergents(a)
    return p[-1], q[-1]

设 \(L_1,L_2,L_3\) 是分式线性变换，并记 \(L_i\) 的系数形成的矩阵为

则它们有如下性质：⁷

1. 分式线性变换的复合 \(L_1\circ L_2\) 和逆变换 \(L_1^{-1}\) 仍然是分式线性变换，即全体分式线性变换构成 群；
2. 分式线性变换在系数同乘以非零常数后保持不变，即对于任意 \(\lambda\neq 0\)，如果 \(M_2=\lambda M_1\)，则 \(L_2=L_1\)；
3. 分式线性变换的复合的系数矩阵，对应着系数矩阵的乘积，即如果 \(M_1M_2=M_3\)，则 \(L_1\circ L_2=L_3\)；
4. 分式线性变换的逆变换的系数矩阵，对应着系数矩阵的逆矩阵，即如果 \(M_1^{-1}=M_2\)，则 \(L_1^{-1}=L_2\)。

此处仅提供分式线性变换的复合和逆变换的形式。得到这个形式后，所有性质都是容易验证的。

分式线性变换 \(L_1\) 和 \(L_2\) 的复合：

分式线性变换 \(L_1(x)\) 的逆变换：

将连分数理解为一列分式线性变换的复合在 \(x=\infty\) 处取值的结果，只需要能够多次查询一段分式线性变换的复合即可。因为每个分式线性变换都可以取逆，所以可以预处理前缀和再用差分的方法查询，复杂度为 \(O(n+m)\) 的；如果需要修改，也可以用线段树等结构存储。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

#include <algorithm>
#include <iostream>
#include <tuple>
#include <vector>

constexpr int M = 1e9 + 7;

// FLTs. Essentially 2x2 matrix.
struct FracLinearTrans {
  int mat[4];

  FracLinearTrans() : mat{} {}

  FracLinearTrans(int x) : mat{x, 1, 1, 0} {}

  FracLinearTrans(int a, int b, int c, int d) : mat{a, b, c, d} {}

  FracLinearTrans operator*(const FracLinearTrans& rhs) const {
    return FracLinearTrans(
        ((long long)mat[0] * rhs.mat[0] + (long long)mat[1] * rhs.mat[2]) % M,
        ((long long)mat[0] * rhs.mat[1] + (long long)mat[1] * rhs.mat[3]) % M,
        ((long long)mat[2] * rhs.mat[0] + (long long)mat[3] * rhs.mat[2]) % M,
        ((long long)mat[2] * rhs.mat[1] + (long long)mat[3] * rhs.mat[3]) % M);
  }

  FracLinearTrans inv() const {
    return FracLinearTrans(mat[3], M - mat[1], M - mat[2], mat[0]);
  }
};

int main() {
  int n, q;
  std::cin >> n >> q;
  // Get prefix sum of FLTs.
  std::vector<FracLinearTrans> ps(1, {1, 0, 0, 1});
  ps.reserve(n + 1);
  for (int i = 1; i <= n; ++i) {
    int a;
    std::cin >> a;
    ps[i] = ps[i - 1] * FracLinearTrans(a);
  }
  // Query.
  for (; q; --q) {
    int l, r;
    std::cin >> l >> r;
    // Difference.
    auto res = ps[l - 1].inv() * ps[r];
    int u = res.mat[0], d = res.mat[2];
    // Correct signs.
    if (!(l & 1)) {
      if (u) u = M - u;
      if (d) d = M - d;
    }
    printf("%d %d\n", u, d);
  }
  return 0;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# PYTHON IS TOO SLOW TO PASS THIS PROBLEM.
# JUST FOR REFERENCE.
M = 10**9 + 7


def mul(a, b):
    return (
        (a[0] * b[0] + a[1] * b[2]) % M,
        (a[0] * b[1] + a[1] * b[3]) % M,
        (a[2] * b[0] + a[3] * b[2]) % M,
        (a[2] * b[1] + a[3] * b[3]) % M,
    )


def inv(a):
    return (a[3], M - a[1], M - a[2], a[0])


n, q = map(int, input().split())
ps = [(1, 0, 0, 1)]
# Get presum.
for a in map(int, input().split()):
    ps.append(mul(ps[-1], (a, 1, 1, 0)))
for _ in range(q):
    l, r = map(int, input().split())
    res = mul(inv(ps[l - 1]), ps[r])
    u, d = res[0], res[2]
    if l % 2 == 0:
        if u:
            u = M - u
        if d:
            d = M - d
    print(u, d)

算法的基本思路就是逐个确定 \(\beta_i\) 的值。记

因为连分数

所以，可以向 \(L\) 通过逐步复合 \(L_{\alpha_k}\) 的方式计算 \(L(\alpha)\) 的大小。但是，如果是希望得到 \(L(\alpha)\) 的连分数表示，那么并不需要完全计算 \(L(\alpha)\) 的值再求出连分数表示。可以在复合 \(L_{\alpha_i}\) 的过程中就能判断 \(\beta_0,\beta_1,\cdots\) 的值。

比如，假设当前计算到

且 \(c_k,d_k\) 同号。那么，\(L\circ L_{\alpha_0}\circ L_{\alpha_1}\circ \cdots \circ L_{\alpha_k}(x)\) 在 \([0,\infty]\) 上单调，且必然取值在 \(\dfrac{a_k}{c_k}\) 和 \(\dfrac{b_k}{d_k}\) 之间。所以，如果

就可以确定它就是 \(L\circ L_{\alpha_0}\circ L_{\alpha_1}\circ \cdots \circ L_{\alpha_k}(x)\) 的整数部分 \(\beta_0\)。此时，在左侧复合 \(L_{\beta_0}^{-1}\) 就可以得到

此时，继续添加 \(L_{\alpha_{k+1}},L_{\alpha_{k+2}},\cdots\) 就可以确定新的整数部分，即 \(\beta_1\)。这样计算下去，直到确定出所有的 \(\beta_j\) 的值。

算法要求 \(c\) 和 \(d\) 同号，是因为要保证函数的不连续点不在 \([0,\infty]\) 范围内。这总是可能的，因为简单连分数的定义要求（除了 \(\alpha_0\) 外的）系数都是正整数。由此可以证明，必然在有限步内成立 \(c\) 和 \(d\) 同号，且将在之后一直保持同号。

具体实现时，只需要维护当前的分式线性变换的系数矩阵 \(\begin{pmatrix}a&b\\c&d\end{pmatrix}\) 并检查 \(c\) 和 \(d\) 是否同号以及 \(\dfrac{a}{c}\) 和 \(\dfrac{b}{d}\) 是否有相同的整数部分。右复合 \(L_{\alpha_i}\) 时，就可以得到 \(\begin{pmatrix}a\alpha_i+b&a\\ c\alpha_i+d&c\end{pmatrix}\)。如果两者整数部分相同为 \(\beta_j\)，则在结果的连分数内添加 \(\beta_j\)，并且左复合 \(L_{\beta_j}^{-1}\)，这相当于计算 \(\begin{pmatrix}c&d\\ a\bmod c & b \bmod d \end{pmatrix}\)。

与单变量的分式线性变换的情形类似，要确定整数部分只需要保证当前的分式线性变换在 \((x,y)\in[0,\infty]\times[0,\infty]\) 内的整数部分保持不变，即 \(e,f,g,h\) 皆同号，且

右复合要替换成计算 \(L(x,y)\mapsto L(L_{\alpha_i}(x),y)\) 和 \(L(x,y)\mapsto L(x,L_{\beta_j}(y))\)，这同样表示成系数的线性变换。左复合则和单变量的情形完全一致，只需要计算取模就可以了。

相较于单变量的情形，双变量的情形需要决定要先复合 \(L_{\alpha_i}\) 还是 \(L_{\beta_j}\)。因为复合的顺序与最后的结果无关，所以可以自由选择复合顺序，比如交替地复合 \(L_{\alpha_i}\) 和 \(L_{\beta_j}\)。或者采用经验法则，优先复合比值差距更大的维度：如果 \(\left|\dfrac{b}{f}-\dfrac{d}{h}\right|>\left|\dfrac{c}{g}-\dfrac{d}{h}\right|\)，那么就先复合 \(L_{\alpha_i}\)；否则，就先复合 \(L_{\beta_j}\)。

循环连分数表示的都是二次无理数。

对于一般的循环连分数 \(x=[a_0,\cdots,a_{k-1},\overline{a_k,\cdots,a_{k+L-1}}]\)，可以设 \(y=[\overline{a_k,\cdots,a_{k+L-1}}]\)，则

其中，\(L_0(\cdot)\) 和 \(L_1(\cdot)\) 都是分式线性变换。于是，得到 \(x\) 满足的方程

不妨设分式线性变换 \(L_0\circ L_1\circ L_0^{-1}(x) = \dfrac{ax+b}{cx+d}\)，则得到 \(x\) 满足的方程

因而，循环连分数都是整系数二次方程的解。又因为无限连分数都是无理数，所以循环连分数都表示了二次无理数。

二次无理数可以表示成循环连分数。

思路是证明余项会重复出现。设二次无理数 \(x\) 可以写作

的形式，其中，\(P_0,Q_0,D\) 都是整数且 \(Q_0\mid D-P_0^2\)。这总是可能的，比如二次无理数 \(x\) 总可以写成

再令 \(P=p_ap_bq_aq_b\)，\(Q=(q_aq_b)^2\) 和 \(D=QD'\) 即可。

将它写成这种形式的好处是，可以证明它的所有余项都具有类似的形式：

其中，\(P_k,Q_k\) 是整数且 \(Q_k\mid D-P_k^2\)。其中，条件 \(Q_k\mid D-P_k^2\) 保证了所有余项的分子中，\(\sqrt{D}\) 前面的系数都是 \(1\)。

为得到余项的形式，可以使用数学归纳法。当 \(k=0\) 时，显然。假设已经得到了 \(r_k\) 的形式，并设 \(a_k=\lfloor r_k\rfloor\)，就有

设 \(r_{k+1}\) 也有类似形式，并和 \(r_k\) 一起代入上式，有

因为二次无理数表示成 \(a+b\sqrt{D}\) 的方式是唯一的，所以比较两侧系数可知

将第二个等式代入第一个等式可以解出 \(P_{k+1}\)：

再代入第二个等式，就可以解出 \(Q_{k+1}\)：

根据归纳假设，\(Q_k\mid D-P_k^2\)，所以确实 \(P_{k+1}\) 和 \(Q_{k+1}\) 都是整数，即 \(r_{k+1}\) 也具有所要求的形式。

最后，证明余项只能取得有限多个值，故而必然重复。前文已经求得余项

而且对于无理数，总有 \(r_k>1\)。同时，它的共轭

对于充分大的 \(k\) 必然小于 \(0\)，因为

这就说明

因此，\(Q_k\) 只能取有限多个值。进而，

所以，\(P_k\) 也只能取有限多个值。故而，余项 \(r_k\) 只有有限多个可能的取值，必然在无限项内重复。

二次无理数总可以表示成

的形式，且 \(Q_0\mid D-P^2_0\)。它的余项

中，\(P_k,Q_k\) 都是整数，且满足递推关系

首先利用上述算法解出 \(\sqrt{x}\) 的周期，将循环节表示成分式线性变换，就可以用 快速幂 获得 \(x_k\) 的值。当然，对于没有进入循环节和不足一个循环节的部分，需要单独处理。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

#include <algorithm>
#include <cmath>
#include <iostream>
#include <tuple>
#include <unordered_map>
#include <vector>

// Return the continued fraction and minimal positive period
//   of a quadratic irrational (x + y * sqrt(n)) / z.
auto quadratic_irrational(int x, int y, int z, int n) {
  int p = x * z;
  int d = n * y * y * z * z;
  int q = z * z;
  int dd = (int)sqrt(n) * y * z;
  int i = 0;
  std::vector<int> a;
  std::unordered_map<size_t, int> used;
  while (!used.count(((1LL << 32) * p) | q)) {
    a.emplace_back((p + dd) / q);
    used[((1LL << 32) * p) | q] = i++;
    p = a.back() * q - p;
    q = (d - p * p) / q;
  }
  return std::make_pair(a, i - used[((1LL << 32) * p) | q]);
}

// Fractional Linear Transformation.
struct FracLinearTrans {
  static constexpr int M = 1e9 + 7;
  int mat[4];

  FracLinearTrans(int a, int b, int c, int d) : mat{a, b, c, d} {}

  FracLinearTrans operator*(const FracLinearTrans& rhs) const {
    return FracLinearTrans(
        ((long long)mat[0] * rhs.mat[0] + (long long)mat[1] * rhs.mat[2]) % M,
        ((long long)mat[0] * rhs.mat[1] + (long long)mat[1] * rhs.mat[3]) % M,
        ((long long)mat[2] * rhs.mat[0] + (long long)mat[3] * rhs.mat[2]) % M,
        ((long long)mat[2] * rhs.mat[1] + (long long)mat[3] * rhs.mat[3]) % M);
  }
};

int main() {
  int x, k, L;
  std::cin >> x >> k;
  std::vector<int> a;
  std::tie(a, L) =
      quadratic_irrational(0, 1, 1, x);  // L==a.size()-1 for sqrt(x)
  FracLinearTrans cyc(1, 0, 0, 1);
  for (int i = a.size() - 1; i; --i) {
    cyc = FracLinearTrans(a[i], 1, 1, 0) * cyc;
  }
  // 1/0=Inf.
  FracLinearTrans res(0, 1, 0, 0);
  // Tail terms.
  for (int i = k % L; i; --i) {
    res = FracLinearTrans(a[i], 1, 1, 0) * res;
  }
  // Binary exponentiation.
  for (int b = k / L; b; b >>= 1) {
    if (b & 1) res = cyc * res;
    cyc = cyc * cyc;
  }
  // First term.
  res = FracLinearTrans(a[0], 1, 1, 0) * res;
  printf("%d/%d", res.mat[1], res.mat[3]);
  return 0;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

from math import sqrt, floor


# Return the continued fraction and minimal positive period
#   of a quadratic irrational (x + y * sqrt(n)) / z.
def quadratic_irrational(x, y, z, n):
    p = x * z
    d = n * y * y * z * z
    q = z * z
    dd = floor(sqrt(n)) * y * z
    i = 0
    a = []
    used = dict()
    while (p, q) not in used:
        a.append((p + dd) // q)
        used[p, q] = i
        p = a[-1] * q - p
        q = (d - p * p) // q
        i += 1
    return a, i - used[p, q]


# Compose (A[0]*x + A[1]) / (A[2]*x + A[3]) and (B[0]*x + B[1]) / (B[2]*x + B[3])
def combine(A, B):
    return [
        t % mod
        for t in [
            A[0] * B[0] + A[1] * B[2],
            A[0] * B[1] + A[1] * B[3],
            A[2] * B[0] + A[3] * B[2],
            A[2] * B[1] + A[3] * B[3],
        ]
    ]


# Binary exponentiation.
def bpow(A, n):
    return (
        [1, 0, 0, 1]
        if not n
        else combine(A, bpow(A, n - 1))
        if n % 2
        else bpow(combine(A, A), n // 2)
    )


mod = 10**9 + 7
x, k = map(int, input().split())
a, T = quadratic_irrational(0, 1, 1, x)

A = (1, 0, 0, 1)  # (x + 0) / (0*x + 1) = x

# apply ak + 1/x = (ak*x+1)/(1x+0) to (Ax + B) / (Cx + D)
for i in reversed(range(1, len(a))):
    A = combine([a[i], 1, 1, 0], A)

C = (0, 1, 0, 0)  # = 1 / 0
while k % T:
    i = k % T
    C = combine([a[i], 1, 1, 0], C)
    k -= 1

C = combine(bpow(A, k // T), C)
C = combine((a[0], 1, 1, 0), C)
print(str(C[1]) + "/" + str(C[3]))

对于纯循环连分数

记

则 \(x\) 和 \(x'\) 互为「倒数负共轭」，即 \(x\) 的共轭的倒数的相反数是 \(x'\)。

因为不要求 \(\ell+1\) 是最小正周期，所以不妨设 \(\ell>0\)。利用反序定理可知，

因为等式两侧都是既约分数，所以

对于纯循环连分数 \(x\)，它的第 \(\ell+1\) 个余项就是 \(x\)，因此，有

所以，它满足二次方程

同理，\(x'\) 满足二次方程

利用系数的关系可知，这个方程可以写作

令 \(y=-\dfrac{1}{x'}\)，则 \(x\) 和 \(y\) 满足同一个方程。但是，\(x>0>y\)，因此它们并非同一个根，而是互为共轭的关系，这就证明了原命题。

二次无理数 \(x\) 可以表示为纯循环连分数，当且仅当 \(x>1\) 且它的共轭 \(-1<x^*<0\)。

如果 \(x\) 是纯循环连分数，那么利用前文的记号，有 \(a_0=a_{\ell+1}\ge 1\)，故而 \(x>1\)。又因为它的倒数负共轭也是循环连分数，所以它的共轭 \(x^*\) 满足 \(-\dfrac{1}{x^*}>1\)，亦即 \(-1<x^*<0\)。这就证明了纯循环连分数都满足该条件。

反过来，设二次无理数 \(x>1\) 且 \(-1<x^*<0\)。对于 \(x\) 的余项 \(r_k\)，有递推关系

等式两边都是二次无理数，取共轭可知

利用这个递推关系，可以证明 \(-1<r_{k}^*<0\) 对所有 \(k\ge 0\) 都成立。

首先，对于 \(k=0\)，有 \(-1<r_0^*=x_0^*<0\)，显然。对于 \(k\ge 0\)，由简单连分数定义和 \(x>1\) 可知，\(a_k\ge 1\)。故而，假设 \(-1<r_k^*<0\)，就有

这就归纳地证明了 \(-1<r_{k}^*<0\) 对所有 \(k\ge 0\) 都成立。因此，有

因为二次无理数一定是循环连分数，所以存在正整数 \(L\)和至少某个充分大的 \(k\)，有 \(r_{k}=r_{k+L}\)。但是，此时必然也有

故而，

这意味着最小的能够使得 \(r_{k}=r_{k+L}\) 成立的 \(k\) 必然是 \(0\)。也就是说，\(x\) 可以表示成纯循环连分数。

对于有理数 \(r>1\)，如果 \(\sqrt{r}\) 是无理数，那么

且对于任意 \(1\le k\le\ell\)，都有 \(a_k = a_{\ell+1-k}\)。

对于二次无理数 \(\sqrt{r}\)，因为 \(\lfloor\sqrt{r}\rfloor+\sqrt{r}>1\) 且 \(-1<\lfloor\sqrt{r}\rfloor-\sqrt{r}<0\)，所以 \(\lfloor\sqrt{r}\rfloor+\sqrt{r}\) 是纯循环连分数：

根据上述定理，它的倒数负共轭具有形式

利用连分数的基本性质可知

但是，又由 \(\lfloor\sqrt{r}\rfloor+\sqrt{r}\) 的连分数表示可知，

因为无理数的连分数表示是唯一的，所以比较中间的系数就知道，\(a_k=a_{\ell+1-k}\) 对所有 \(1\le k\le\ell\) 都成立。

对于无界集合 \(x\ge 0\)，上凸壳就是直线 \(y=rx\) 本身。然而，如下图所示，如果还要求 \(x\le N\)，那么上凸壳最终会偏离直线。

从 \((0,0)\) 开始，可以自左向右地求出上凸壳的所有整点。假设当前已经求出的上凸壳的最后一个整点是 \((x,y)\)。现在要求出下一个整点 \((x',y')\)。顶点 \((x',y')\) 在 \((x,y)\) 右上方，记 \((\Delta x,\Delta y)=(x'-x,y'-y)\) 为两者的差值。那么，必然有

第二个不等式成立，因为条件 \(\Delta y>r\Delta x\) 与 \((x,y)\) 已经在上凸壳上这件事矛盾。观察 \((\Delta x,\Delta y)\) 需要满足的条件，对于不同的点 \((x,y)\)，只有 \(\Delta x\) 的上界在变化。所以，只要能解决这个子问题，就可以递归地求出原问题的所有整点。

进而，考虑子问题的解法。对比于原问题，子问题相当于将 \(x\) 的上界修改为 \(N'\)，并求出上凸壳中与原点相邻的第一个整点。记子问题的解为 \((q,p)\)。那么，\(p\) 与 \(q\) 必然是互素的（否则不是第一个整点），且与原点连线的斜率 \(\dfrac{p}{q}\) 是所有位于直线 \(y=rx\) 下方且横坐标不超过 \(N'\) 的整点中最大的（否则不在凸包上）。结合前文的 几何解释 可知，这样的点 \((x,y)\) 必然对应于 \(r\) 的一个下中间分数。因为分母越大的下中间分数离 \(r\) 越近，所以子问题的解 \((q,p)\) 对应着所有分母不超过 \(N'\) 的下中间分数中分母最大的那个。

当然，实际求解时，没必要对每个子问题都重新求出这样的下中间分数。应该首先求出所有的渐近分数，这相当于提供了遍历所有的下中间分数的方法。然后分母从大到小地遍历下中间分数，每次都尝试将它加到前一个整点 \((x,y)\) 上，直到不能添加为止才继续尝试下一个下中间分数。

此处有一些显然的优化。首先，对于下中间分数 \((q,p)\)，必然存在奇数 \(k\) 和 \(0\le t<a_k\) 使得 \((q,p)=(q_{k-1},p_{k-1})+t(q_k,p_k)\)。只要找到最大的 \(t\) 使得 \(q_{k-1}+tq_k+x\le N\) 满足就好了，亦即 \(t=\left\lfloor\dfrac{N-q_{k-1}-x}{q_k}\right\rfloor\)。不用担心 \(t\) 越界，因为更大的下渐近分数 \((q_{k+2},p_{k+2})\) 已经添加完了。而每次确定添加的次数的时候，直接计算 \(\left\lfloor\dfrac{N-x}{q}\right\rfloor\) 即可，不必逐个尝试。

优化后的算法的复杂度是 \(O(n)\) 的。虽然下中间分数对应的整点可能有很多，但是真正成为增量的并不多。下面要说明，所有 \(0\le t<a_k\) 的下中间分数 \((q,p)=(q_{k-1},p_{k-1})+t(q_k,p_k)\) 中，至多会出现两个增量。假设这些下中间分数中确实出现了增量，则此时必然有 \(q_{k-1}\le N-x<q_{k+1}\)。不妨设 \(t=\left\lfloor\dfrac{N-q_{k-1}-x}{q_k}\right\rfloor\)。如果 \(t=0\)，则增量就有 \(\Delta x=q_{k-1}\)，故而添加完增量后，就有 \(N-x'<q_{k-1}\)，不会再在这些下中间分数中出现新的增量；如果 \(t>0\)，那么添加完增量后，必然有 \(N-x'=(N-q_{k-1}-x)\bmod q_k<q_k\)，即使还会在同一段下中间分数中出现新的增量，下次也只能有 \(t'=0\)。因此，在这样的一段下中间分数中，至多只能出现两个增量。这就说明，总的时间复杂度是 \(O(n)\) 的。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

// Find [ah, ph, qh] such that points r[i]=(ph[i], qh[i]) constitute
// upper convex hull of lattice points on 0 <= x <= N and 0 <= y <= r * x,
// where r = [a0, a1, a2, ...] and there are ah[i]-1 integer points on the
// segment between r[i] and r[i+1].
auto hull(std::vector<int> a, int N) {
  std::vector<int> p, q;
  std::tie(p, q) = convergents(a);
  int t = N / q.back();
  std::vector<int> ah = {t};
  std::vector<int> ph = {0, t * p.back()};
  std::vector<int> qh = {0, t * q.back()};

  for (int i = q.size() - 1; i; --i) {
    if (i % 2) {
      while (qh.back() + q[i - 1] <= N) {
        t = (N - qh.back() - q[i - 1]) / q[i];
        int dp = p[i - 1] + t * p[i];
        int dq = q[i - 1] + t * q[i];
        int k = (N - qh.back()) / dq;
        ah.push_back(k);
        ph.push_back(ph.back() + k * dp);
        qh.push_back(qh.back() + k * dq);
      }
    }
  }
  return make_tuple(ah, ph, qh);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

# Find [ah, ph, qh] such that points r[i]=(ph[i], qh[i]) constitute
# upper convex hull of lattice points on 0 <= x <= N and 0 <= y <= r * x,
# where r = [a0, a1, a2, ...] and there are ah[i]-1 integer points on the
# segment between r[i] and r[i+1].
def hull(a, N):
    p, q = convergents(a)
    t = N // q[-1]
    ah = [t]
    ph = [0, t * p[-1]]
    qh = [0, t * q[-1]]
    for i in reversed(range(len(q))):
        if i % 2 == 1:
            while qh[-1] + q[i - 1] <= N:
                t = (N - qh[-1] - q[i - 1]) // q[i]
                dp = p[i - 1] + t * p[i]
                dq = q[i - 1] + t * q[i]
                k = (N - qh[-1]) // dq
                ah.append(k)
                ph.append(ph[-1] + k * dp)
                qh.append(qh[-1] + k * dq)
    return ah, ph, qh

这个问题有一个复杂度为 \(O(\sqrt N)\) 的解法：不妨设 \(A\ge B\)，因为 \(A(B+x)+By=Ax+B(A+y)\)，所以只需要在 \(x\le\min\{N/A, B\}\) 中搜索答案即可。这足够通过本题。但是，如果应用连分数方法，那么时间复杂度就可以降低到 \(O(\log N)\)。

为了讨论方便，首先通过代换 \(x\mapsto\left\lfloor N/A\right\rfloor-x\) 来改变 \(x\) 的符号。令 \(C=N\bmod A\) 和 \(M=\left\lfloor N/A\right\rfloor\)，则原问题转化为在 \(0\le x\le M\) 且 \(By-Ax\le C\) 的条件下，求最优的 \((x,y)\) 使得 \(By-Ax\) 最大。对于每个固定的 \(x\)，最优的 \(y\) 的取值为 \(\left\lfloor\dfrac{Ax+C}{B}\right\rfloor\)。

接下来要说明的是，这个问题和上一个例题具有类似的解法。但是，与上一个例题中使用下中间分数偏离直线不同，本题需要使用上中间分数来接近直线。具体来说，\(C-(By-Ax)\) 的值正比于点 \((x,y)\) 与直线 \(By-Ax=C\) 的距离。要最大化 \(By-Ax\)，就等价于最小化这个距离。算法的目标是要找到直线 \(By-Ax=C\) 下方距离它最近的可行的整点。算法的思路就是从最左侧的点开始，沿着这些整点的上凸壳搜索，逐步缩小与直线的距离，直到得到最优解。

在 \((x,y)\) 的坐标系内，算法从 \((0,\lfloor C/B\rfloor)\) 出发，递归地寻找并添加最优的增量 \((\Delta x,\Delta y)\)，且保证添加后的点比起之前更靠近直线 \(By-Ax=C\)，但是不能到达直线的另一侧，也不能让横坐标大于 \(M\)。设已经得到的点是 \((x,y)\)，那么增量 \((\Delta x,\Delta y)\) 满足的条件就是

按照沿下凸壳搜索的思路，只需要找到满足这些条件的点中 \(\Delta x\) 最小的即可。将第一个不等式改写成

结合前文的 几何解释 可知，只要后面的常数项小于 \(1\)，那么满足这个不等式的整点 \((\Delta x,\Delta y)\) 中横坐标最小的，一定对应着某个上中间分数。这是因为它是所有分母不超过它的分母的分数中，从上方逼近某个实数效果最好的，这只能是上中间分数。而每次添加增量后，都会导致 \(\Delta y\) 的上界变得更紧，这意味着必须考察分母更大的上中间分数。

仿照上一个例题的思路。分母从小到大考察所有上中间分数，如果能够找到横坐标和纵坐标都不越界的上中间分数，就添加进去，并更新相应的上界。当所有可行的上中间分数都添加结束后，得到的就是最优解。相较于之前，这个题目需要同时保证横纵坐标都不越界，需要格外注意。基于和上一个例题类似的论述，不过这次是使用 \(B\Delta y-A\Delta x\) 代替之前的 \(\Delta x\)，可以说明这个算法的复杂度是 \(O(\log\min\{A,B\})\) 的。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

// Find (x, y) such that y = (A*x+B)/C,
// such that Cy - Ax is max and 0 <= x <= N.
auto closest(int A, int B, int C, int N) {
  // y <= (A*x + B)/C <=> diff(x, y) <= B
  auto diff = [&](int x, int y) -> int { return C * y - A * x; };
  std::vector<int> p, q;
  std::tie(p, q) = convergents(fraction(A, C));
  int qh = 0, ph = B / C;
  for (int i = 2; i < q.size() - 1; ++i) {
    if (i % 2 == 0) {
      while (diff(qh + q[i + 1], ph + p[i + 1]) <= B) {
        int t = 1 + (diff(qh + q[i - 1], ph + p[i - 1]) - B - 1) /
                        (-diff(q[i], p[i]));
        int dp = p[i - 1] + t * p[i];
        int dq = q[i - 1] + t * q[i];
        int k = (N - qh) / dq;
        if (k == 0) {
          return std::make_pair(qh, ph);
        }
        if (diff(dq, dp) != 0) {
          k = std::min(k, (B - diff(qh, ph)) / diff(dq, dp));
        }
        qh += k * dq;
        ph += k * dp;
      }
    }
  }
  return std::make_pair(qh, ph);
}

auto solve(int A, int B, int N) {
  int x, y;
  std::tie(x, y) = closest(A, N % A, B, N / A);
  return std::make_pair(N / A - x, y);
}