类欧几里德算法

引入

类欧几里德算法是洪华敦在 2016 年冬令营营员交流中提出的内容。它常用于解决形如

⌊ \frac{a i + b}{c} ⌋

结构的数列（下标为 $i$ ）的求和问题。它的主要想法是，利用分数自身的递归结构，将问题转化为更小规模的问题，递归求解。因为分数的递归结构和欧几里得算法存在直接的联系，因此，这一求和方法也称为类欧几里得算法。

因为连分数和 Stern–Brocot 树等方法同样刻画了分数的递归结构，所以利用类欧几里得算法可以解决的问题，通常也可以用这些方法解决。与这些方法相比，类欧几里得算法通常更容易理解，它的实现也更为简明。

类欧几里得算法

最简单的例子，就是求和问题：

f (a, b, c, n) = \sum_{i = 0}^{n} ⌊ \frac{a i + b}{c} ⌋,

其中， $a, b, c, n$ 都是正整数。

代数解法

首先，将 $a, b$ 对 $c$ 取模，可以简化问题，将问题转化为 $0 \leq a, b < c$ 的情形：

\begin{aligned} f (a, b, c, n) & = \sum_{i = 0}^{n} ⌊ \frac{a i + b}{c} ⌋ \\ = \sum_{i = 0}^{n} ⌊ \frac{(⌊ \frac{a}{c} ⌋ c + (a mod c)) i + (⌊ \frac{b}{c} ⌋ c + (b mod c))}{c} ⌋ \\ = \sum_{i = 0}^{n} (⌊ \frac{a}{c} ⌋ i + ⌊ \frac{b}{c} ⌋ + ⌊ \frac{(a mod c) i + (b mod c)}{c} ⌋) \\ = \frac{n (n + 1)}{2} ⌊ \frac{a}{c} ⌋ + (n + 1) ⌊ \frac{b}{c} ⌋ + f (a mod c, b mod c, c, n) . \end{aligned}

现在，考虑转化后的问题。令

m = ⌊ \frac{a n + b}{c} ⌋ .

那么，原问题可以写作二次求和式：

\sum_{i = 0}^{n} ⌊ \frac{a i + b}{c} ⌋ = \sum_{i = 0}^{n} \sum_{j = 0}^{m - 1} [j < ⌊ \frac{a i + b}{c} ⌋] .

交换求和次序，这需要对于每个 $j$ 计算满足条件的 $i$ 的范围。为此，将条件变形：

\begin{aligned} j < ⌊ \frac{a i + b}{c} ⌋ = ⌈ \frac{a i + b + 1}{c} ⌉ - 1 \\ ⟺ j + 1 < ⌈ \frac{a i + b + 1}{c} ⌉ ⟺ j + 1 < \frac{a i + b + 1}{c} \\ ⟺ \frac{c j + c - b - 1}{a} < i ⟺ ⌊ \frac{c j + c - b - 1}{a} ⌋ < i . \end{aligned}

变形过程中多次利用了取整函数的性质。代入变形后的条件，原式可以写作：

\begin{aligned} f (a, b, c, n) & = \sum_{j = 0}^{m - 1} \sum_{i = 0}^{n} [i > ⌊ \frac{c j + c - b - 1}{a} ⌋] \\ = \sum_{j = 0}^{m - 1} (n - ⌊ \frac{c j + c - b - 1}{a} ⌋) \\ = n m - f (c, c - b - 1, a, m - 1) . \end{aligned}

令 $(a^{'}, b^{'}, c^{'}, n^{'}) = (c, c - b - 1, a, m - 1)$ ，这就又回到了前面讨论过的 $a^{'} > c^{'}$ 的情形。

将这两步转化结合在一起，可以发现在过程中， $(a, c)$ 不断地取模后交换位置，直到 $a = 0$ 。这就类似于对 $(a, c)$ 进行辗转相除，这也是类欧几里德算法的得名。它的时间复杂度是 $O (\log min {a, c})$ 的。

在计算过程中，可能会出现 $m = 0$ 的情形，此时内层递归会出现 $n = - 1$ 。这并不影响最终的结果。但是，如果要求出现 $m = 0$ 时，直接终止算法，那么算法的时间复杂度可以改良为 $O (\log min {a, c, n})$ 的。

对复杂度的解释

利用该算法和欧几里得算法的相似性，很容易说明它的时间复杂度是 $O (\log min {a, c})$ 的。因此，只需要说明，如果在 $m = 0$ 时终止算法，那么它的时间复杂度也是 $O (\log n)$ 的。

令 $m = ⌊ (a n + b) / c ⌋$ ，并记 $S = m n$ ， $k = m / n$ ，它们分别相当于几何直观（见下一节）中点阵图的面积和直线的斜率。对于充分大的 $n$ ，近似有 $k ≐ a / c$ 。

考察 $S$ 和 $k$ 在算法过程中的变化。第一步取模时， $n$ 保持不变， $k$ 近似由 $a / c$ 变为 $(a mod c) / c$ ，相当于斜率由 $k$ 变为 $k - ⌊ k ⌋$ ，而 $S$ 也近似变为原来的 $(k - ⌊ k ⌋)$ 倍。第二步交换横纵坐标时， $S$ 近似保持不变， $k$ 则变为它的倒数。因此，若设两步操作后，二元组 $(k, S)$ 变为 $(k^{'}, S^{'})$ ，则有 $k^{'} = (k - ⌊ k ⌋)^{- 1}$ 且 $S^{'} = (k - ⌊ k ⌋) S$ 。

因为 $1 \leq ⌊ k^{'} ⌋ \leq k^{'} < ⌊ k^{'} ⌋ + 1$ ，所以，递归计算两轮后，乘积缩小的倍数最少为

(k^{'} - ⌊ k^{'} ⌋) (k - ⌊ k ⌋) = 1 - \frac{⌊ k^{'} ⌋}{k^{'}} < 1 - \frac{⌊ k^{'} ⌋}{⌊ k^{'} ⌋ + 1} = \frac{1}{⌊ k^{'} ⌋ + 1} \leq \frac{1}{2} .

因此，至多 $O (\log S)$ 轮，算法必然终止。因为从第二轮开始，每轮开始时的 $S$ 总是不超过上一轮取模结束后的 $S$ ，而后者大致为 $k n^{2}$ 且 $k < 1$ ，因而 $O (\log S) \subseteq O (\log n)$ 。这就得到了上述结论。

模板题的参考实现如下：

模板题实现（Library Checker - Sum of Floor of Linear）

#include <iostream>

long long solve(long long a, long long b, long long c, long long n) {
  long long n2 = n * (n + 1) / 2;
  if (a >= c || b >= c)
    return solve(a % c, b % c, c, n) + (a / c) * n2 + (b / c) * (n + 1);
  long long m = (a * n + b) / c;
  if (!m) return 0;
  return m * n - solve(c, c - b - 1, a, m - 1);
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    int a, b, c, n;
    std::cin >> n >> c >> a >> b;
    std::cout << solve(a, b, c, n - 1) << '\n';
  }
  return 0;
}

几何直观

这个算法还可以从几何的角度理解。类欧几里得算法可以解决的问题主要是直线下整点计数问题。

如下图最左部分所示，该求和式相当于求直线

y = \frac{a x + b}{c}

下方， $x$ 轴上方（不包括 $x$ 轴），且横坐标位于 $[0, n]$ 之间的格点数目。

首先，移除斜率和截距中的整数部分。这一步相当于将上图中间部分的蓝点数量单独计算出来。当斜率和截距都是整数时，蓝点一定构成一个梯形阵列，也就是说，不同纵列的格点形成等差数列，因而这些点的数量是容易计算的。将这些点移除后，剩余的格点和上图最右部分的红点数量一致。问题就转化成了斜率和截距都小于一的情形。因为梯形的高为 $n + 1$ 且两个底边长度分别为 $⌊ b / c ⌋$ 和 $(⌊ a / c ⌋ n + ⌊ b / c ⌋)$ ，所以，利用梯形面积公式，这一步骤可以归纳为算式

f (a, b, c, n) = f (a mod c, b mod c, c, n) + \frac{1}{2} (n + 1) (⌊ \frac{b}{c} ⌋ + (⌊ \frac{a}{c} ⌋ n + ⌊ \frac{b}{c} ⌋)) .

然后，翻转横纵坐标轴。如下图最左部分所示，图中的红点和蓝点构成了一个横向长度为 $n$ 、纵向长度为 $m = ⌊ (a n + b) / c ⌋$ 的矩形点阵。要计算红点的数量，只需要计算蓝点的数量，再用矩形点阵的数量减去蓝点的数量即可。翻转后，上图左半部分中的蓝点点阵就变成了某条直线下的红色点阵。而且，翻转后，斜率大于一，就又回到了上文已经处理过的情形。

关键在于如何计算新的红色点阵上方的直线的方程。将上图最左部分的横纵坐标轴翻转，就得到上图中间部分。翻转后的红色点阵上方的直线（中间部分的实线），并非对应翻转前的直线（最左部分的实线），而是翻转前的直线向左上平移一点点的结果（最左部分的虚线）。这是因为，如果直接将直线（最左部分的实线）翻转，将得到中间部分的虚线，但是按照定义，它下方的格点包括恰好落在直线上的格点，这就会导致直线上的格点重复计数。为了避免这一点，需要将翻转直线 $y = (a x + b) / c$ 后得到的直线 $y = (c x - b) / a$ 向下平移一点点，得到直线 $y = (c x - b - 1) / a$ ，这样它下方的点阵才恰为翻转前的蓝色点阵。

还有另一处细节需要处理。上图中间部分的直线的截距是负数，这意味着还没有回到之前的初始情形。要让截距恢复为非负数，只需要将直线（中间部分的实线）向左平移一个单位。这样做不会漏掉任何格点，因为翻转前的蓝色点阵中没有纵坐标为零的点，翻转后也就不存在横坐标为零的点。最后，直线方程就变为 $y = (c x + c - b - 1) / a$ ；同时，点阵的横坐标的上界也从 $m$ 变成了 $m - 1$ 。这一步骤可以归纳为算式

f (a, b, c, n) = m n - f (c, c - b - 1, a, m - 1) .

这种递归的算法行得通，主要有两个原因：

第一，直线的斜率不断地先取小数部分再取倒数，这等价于计算直线斜率 $k = a / c$ 的连分数展开。因为有理分数的连分数展开的长度是 $O (\log min {a, c})$ 的，所以这一过程一定在 $O (\log min {a, c})$ 步后终止；
第二，因为每次翻转坐标轴的时候，直线斜率都是小于一的，因此，直觉上应该有 $m < n$ ，也就是说，经过这样一轮迭代后，横坐标的范围一直是在缩小的。前文的复杂度计算中通过严格的分析说明，每两轮迭代后， $n$ 至多为原来的一半，故而这一过程一定在 $O (\log n)$ 步后终止。

这也是斜率为有理数时的类欧几里得算法的复杂度是 $O (\log min {a, c, n})$ 的原因。

利用类似的几何直观，可以将类欧几里得算法推广到斜率为无理数的情形，具体分析请参考后文的例题。

例题

【模板】类欧几里得算法

多组询问。给定正整数 $a, b, c, n$ ，求

\begin{aligned} f (a, b, c, n) & = \sum_{i = 0}^{n} ⌊ \frac{a i + b}{c} ⌋, \\ g (a, b, c, n) & = \sum_{i = 0}^{n} i ⌊ \frac{a i + b}{c} ⌋, \\ h (a, b, c, n) & = \sum_{i = 0}^{n} {⌊ \frac{a i + b}{c} ⌋}^{2} . \end{aligned}

解答一

类似于 $f$ 的推导，可以得到 $g, h$ 的递归表达式。

首先，利用取模，将问题转化为 $0 \leq a, b < c$ 的情形：

\begin{aligned} g (a, b, c, n) & = g (a mod c, b mod c, c, n) + ⌊ \frac{a}{c} ⌋ \frac{n (n + 1) (2 n + 1)}{6} + ⌊ \frac{b}{c} ⌋ \frac{n (n + 1)}{2}, \\ h (a, b, c, n) & = h (a mod c, b mod c, c, n) \\ + 2 ⌊ \frac{b}{c} ⌋ f (a mod c, b mod c, c, n) + 2 ⌊ \frac{a}{c} ⌋ g (a mod c, b mod c, c, n) \\ + {⌊ \frac{a}{c} ⌋}^{2} \frac{n (n + 1) (2 n + 1)}{6} + {⌊ \frac{b}{c} ⌋}^{2} (n + 1) + ⌊ \frac{a}{c} ⌋ ⌊ \frac{b}{c} ⌋ n (n + 1) . \end{aligned}

然后，利用交换求和次序，可以进一步转化。同样地，令

m = ⌊ \frac{a n + b}{c} ⌋ .

那么，对于和式 $g$ ，有

\begin{aligned} g (a, b, c, n) & = \sum_{i = 0}^{n} i ⌊ \frac{a i + b}{c} ⌋ \\ = \sum_{i = 0}^{n} \sum_{j = 0}^{m - 1} i [j < ⌊ \frac{a i + b}{c} ⌋] \\ = \sum_{j = 0}^{m - 1} \sum_{i = 0}^{n} i [i > ⌊ \frac{c j + c - b - 1}{a} ⌋] \\ = \sum_{j = 0}^{m - 1} \frac{1}{2} (⌊ \frac{c j + c - b - 1}{a} ⌋ + n + 1) (n - ⌊ \frac{c j + c - b - 1}{a} ⌋) \\ = \frac{1}{2} m n (n + 1) - \frac{1}{2} \sum_{j = 0}^{m - 1} ⌊ \frac{c j + c - b - 1}{a} ⌋ - \frac{1}{2} \sum_{j = 0}^{m - 1} {⌊ \frac{c j + c - b - 1}{a} ⌋}^{2} \\ = \frac{1}{2} m n (n + 1) - \frac{1}{2} f (c, c - b - 1, a, m - 1) - \frac{1}{2} h (c, c - b - 1, a, m - 1) . \end{aligned}

对于和式 $h$ ，有

\begin{aligned} h (a, b, c, n) & = \sum_{i = 0}^{n} {⌊ \frac{a i + b}{c} ⌋}^{2} \\ = \sum_{i = 0}^{n} \sum_{j = 0}^{m - 1} (2 j + 1) [j < ⌊ \frac{a i + b}{c} ⌋] \\ = \sum_{j = 0}^{m - 1} \sum_{i = 0}^{n} (2 j + 1) [i > ⌊ \frac{c j + c - b - 1}{a} ⌋] \\ = \sum_{j = 0}^{m - 1} (2 j + 1) (n - ⌊ \frac{c j + c - b - 1}{a} ⌋) \\ = n m^{2} - \sum_{j = 0}^{m - 1} ⌊ \frac{c j + c - b - 1}{a} ⌋ - 2 \sum_{j = 0}^{m - 1} j ⌊ \frac{c j + c - b - 1}{a} ⌋ \\ = n m^{2} - f (c, c - b - 1, a, m - 1) - 2 g (c, c - b - 1, a, m - 1) . \end{aligned}

从几何直观的角度看，这些非线性的求和式相当于给区域中的每个点 $(i, j)$ 都赋予了相应的权重 $w (i, j)$ 。除了这些权重之外，其余部分的计算过程是完全一致的。对于权重的选择，一般地，有

\sum_{i = 0}^{n} i^{r} {⌊ \frac{a i + b}{c} ⌋}^{s} = \sum_{i = 0}^{n} \sum_{j = 0}^{m - 1} i^{r} ((j + 1)^{s} - j^{s}) [j < ⌊ \frac{a i + b}{c} ⌋] .

本题的另一个特点是， $g$ 和 $h$ 在递归计算时，会相互交错。因此，需要将 $(f, g, h)$ 作为三元组同时递归。

#include <iostream>

struct Data {
  int f, g, h;
};

Data solve(long long a, long long b, long long c, long long n) {
  constexpr long long M = 998244353;
  constexpr long long i2 = (M + 1) / 2;
  constexpr long long i6 = (M + 1) / 6;
  long long n2 = (n + 1) * n % M * i2 % M;
  long long n3 = (2 * n + 1) * (n + 1) % M * n % M * i6 % M;
  Data res = {0, 0, 0};
  if (a >= c || b >= c) {
    auto tmp = solve(a % c, b % c, c, n);
    long long aa = a / c, bb = b / c;
    res.f = (tmp.f + aa * n2 + bb * (n + 1)) % M;
    res.g = (tmp.g + aa * n3 + bb * n2) % M;
    res.h = (tmp.h + 2 * bb * tmp.f % M + 2 * aa * tmp.g % M +
             aa * aa % M * n3 % M + bb * bb % M * (n + 1) % M +
             2 * aa * bb % M * n2 % M) %
            M;
    return res;
  }
  long long m = (a * n + b) / c;
  if (!m) return res;
  auto tmp = solve(c, c - b - 1, a, m - 1);
  res.f = (m * n - tmp.f + M) % M;
  res.g = (m * n2 + (M - tmp.f) * i2 + (M - tmp.h) * i2) % M;
  res.h = (n * m % M * m - tmp.f - tmp.g * 2 + 3 * M) % M;
  return res;
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    int n, a, b, c;
    std::cin >> n >> a >> b >> c;
    auto res = solve(a, b, c, n);
    std::cout << res.f << ' ' << res.h << ' ' << res.g << '\n';
  }
  return 0;
}

[清华集训 2014] Sum

多组询问。给定正整数 $n$ 和 $r$ ，求

\sum_{d = 1}^{n} (- 1)^{⌊ d \sqrt{r} ⌋} .

解答一

如果 $r$ 是完全平方数，那么当 $\sqrt{r}$ 为偶数时，和式为 $n$ ；否则，和式依据 $n$ 奇偶性不同，在 $0$ 和 $- 1$ 之间交替变化。下面考虑 $r$ 不是完全平方数的情形。

为了应用类欧几里得算法，首先将求和式转化为熟悉的形式：

\begin{aligned} \sum_{d = 1}^{n} (- 1)^{⌊ d \sqrt{r} ⌋} & = \sum_{d = 1}^{n} (1 - 2 (⌊ d \sqrt{r} ⌋ mod 2)) \\ = n - 2 \sum_{d = 1}^{n} (⌊ d \sqrt{r} ⌋ - 2 ⌊ \frac{⌊ d \sqrt{r} ⌋}{2} ⌋) \\ = n - 2 \sum_{d = 1}^{n} ⌊ d \sqrt{r} ⌋ + 4 \sum_{d = 1}^{n} ⌊ \frac{d \sqrt{r}}{2} ⌋ \\ = n - 2 f (n, 1, 0, 1) + 4 f (n, 1, 0, 2) \end{aligned}

其中的函数 $f$ 具有形式

f (a, b, c, n) = \sum_{i = 1}^{n} ⌊ \frac{a \sqrt{r} + b}{c} i ⌋ .

与正文中的算法不同的是，此处的斜率不再是有理数。设斜率

k = \frac{a \sqrt{r} + b}{c} .

同样分为两种情形讨论。如果 $k \geq 1$ ，那么

\begin{aligned} f (a, b, c, n) & = \sum_{i = 1}^{n} ⌊ k i ⌋ = \sum_{i = 1}^{n} ⌊ (k - ⌊ k ⌋) i ⌋ + ⌊ k ⌋ \sum_{i = 1}^{n} i \\ = ⌊ k ⌋ \frac{n (n + 1)}{2} + f (a, b - c ⌊ k ⌋, c, n) . \end{aligned}

问题转化为斜率小于一的情形。如果 $k < 1$ ，那么设 $m = ⌊ n k ⌋$ ，有

\begin{aligned} f (a, b, c, n) & = \sum_{i = 1}^{n} ⌊ k i ⌋ = \sum_{i = 1}^{n} \sum_{j = 1}^{m} [j \leq ⌊ k i ⌋] \\ = \sum_{j = 1}^{m} \sum_{i = 1}^{n} [i > ⌊ k^{- 1} j ⌋] = n m - \sum_{j = 1}^{m} \sum_{i = 1}^{n} [i \leq ⌊ k^{- 1} j ⌋] . \end{aligned}

此处的推导中，交换 $i$ 和 $j$ 的条件比正文中的情形更为简单，是因为直线 $y = k x$ 上没有除了原点之外的格点。关键在于交换后的求和式写成 $f (a, b, c, n)$ 的形式，这相当于要求 $a^{'}, b^{'}, c^{'}$ 满足

k^{- 1} = \frac{a^{'} \sqrt{r} + b^{'}}{c^{'}} .

这并不困难，只需要将分母有理化，就能得到

k^{- 1} = \frac{c}{a \sqrt{r} + b} = \frac{c a \sqrt{r} - c b}{a^{2} r - b^{2}} .

因此，有

a^{'} = c a, b^{'} = - c b, c^{'} = a^{2} r - b^{2} .

这说明

f (a, b, c, n) = n m - f (c a, - c b, a^{2} r - b^{2}, m) .

为了避免整数溢出，需要每次都将 $a, b, c$ 同除以它们的最大公约数。因为这个计算过程和计算 $k$ 的连分数的过程完全一致，所以根据连分数理论，只要保证 $gcd (a, b, c) = 1$ ，它们在计算过程中必然在整型范围内。另外，尽管 $(a, b, c, n)$ 不会溢出，但是在该题数据范围下， $f (a, b, c, n)$ 可能会超过 $64$ 位整数的范围，自然溢出即可，无需额外处理，最后结果一定在 $[- n, n]$ 之间。

尽管斜率不会变为零，算法的复杂度仍然是 $O (\log n)$ 的，这一点从前文关于算法复杂度的论证容易看出。

#include <cmath>
#include <iostream>

long long r;
long double sqrt_r;

long long gcd(long long a, long long b) { return b ? gcd(b, a % b) : a; }

unsigned long long f(long long a, long long b, long long c, long long n) {
  if (!n) return 0;
  auto d = gcd(a, gcd(b, c));
  a /= d;
  b /= d;
  c /= d;
  unsigned long long k = (a * sqrt_r + b) / c;
  if (k) {
    return n * (n + 1) / 2 * k + f(a, b - c * k, c, n);
  } else {
    unsigned long long m = n * (a * sqrt_r + b) / c;
    return n * m - f(c * a, -c * b, a * a * r - b * b, m);
  }
}

unsigned long long solve(long long n, long long r) {
  long long sqr = sqrt_r = sqrtl(r);
  if (r == sqr * sqr) return r % 2 ? (n % 2 ? -1 : 0) : n;
  return n - 2 * f(1, 0, 1, n) + 4 * f(1, 0, 2, n);
}

int main() {
  int t;
  std::cin >> t;
  for (; t; --t) {
    int n;
    std::cin >> n >> r;
    long long res = solve(n, r);
    std::cout << res << '\n';
  }
  return 0;
}

Fraction

给定正整数 $a, b, c, d$ ，求所有满足 $a / b < p / q < c / d$ 的最简分数 $p / q$ 中 $(q, p)$ 的字典序最小的那个。

解答

这道题目也是 Stern–Brocot 树的经典应用，相关题解可以在此处找到。因为它只依赖于分数的递归结构，所以它同样可以利用类似欧几里得算法的方法求解，故而也可以视作类欧几里得算法的一个应用。

如果 $a / b$ 和 $c / d$ 之间（不含端点）存在至少一个自然数，可以直接取 $(q, p) = (1, ⌊ a / b ⌋ + 1)$ 。否则，必然有

⌊ \frac{a}{b} ⌋ \leq \frac{a}{b} < \frac{p}{q} < \frac{c}{d} \leq ⌊ \frac{a}{b} ⌋ + 1.

从这个不等式中可以看出， $p / q$ 的整数部分可以确定为 $⌊ a / b ⌋$ ，直接消去该整数部分，然后整体取倒数，用于确定它的小数部分。这正是确定 $p / q$ 的连分数的基本方法。若最终的答案是 $p / q$ ，那么算法的时间复杂度为 $O (\log min {p, q})$ 。

此处，有一个细节问题，即取倒数之后得到的字典序最小的分数，是否是取倒数之前的字典序最小的分数。换句话说，满足 $a / b < p / q < c / d$ 的分数 $p / q$ 中，字典序 $(q, p)$ 最小的，是否也是字典序 $(p, q)$ 最小的。假设不然，设 $p / q$ 是字典序 $(q, p)$ 最小的，但是 $r / s \neq p / q$ 是字典序 $(r, s)$ 最小的。这必然有 $r < p$ 且 $q < s$ 。但是，这说明

\frac{a}{b} < \frac{r}{s} < \frac{r}{q} < \frac{p}{q} < \frac{c}{d} .

因此， $r / q$ 无论按照哪个字典序怎样都是严格更小于当前解的。这与所设条件矛盾。因此，上述算法是正确的。

#include <iostream>

void solve(int a, int b, int& p, int& q, int c, int d) {
  if ((a / b + 1) * d < c) {
    p = a / b + 1;
    q = 1;
  } else {
    solve(d, c - d * (a / b), q, p, b, a % b);
    p += q * (a / b);
  }
}

int main() {
  int a, b, c, d, p, q;
  while (std::cin >> a >> b >> c >> d) {
    solve(a, b, p, q, c, d);
    std::cout << p << '/' << q << '\n';
  }
  return 0;
}

万能欧几里得算法

上一节讨论的类欧几里得算法推导通常较为繁琐，而且能够解决的和式主要是可以转化为直线下（带权）整点计数问题的和式。本节讨论一种更为一般的方法，它进一步抽象了上述过程，从而可以解决更多的问题。因此，这一方法也称为万能欧几里得算法。它同样利用了分数的递归结构求解问题，但是与类欧几里得算法约化问题的思路稍有不同。

仍然考虑最经典的求和问题：

f (a, b, c, n) = \sum_{i = 1}^{n} ⌊ \frac{a i + b}{c} ⌋,

其中， $a, b, c, n$ 都是正整数。

问题转化

设参数为 $(a, b, c, n)$ 的线段为

y = \frac{a x + b}{c}, 0 < x \leq n .

对于这条线段，可以按照如下方法定义一个由 $U$ 和 $R$ 组成的字符串 $S$ ，也称为 操作序列：

字符串恰有 $n$ 个 $R$ 和 $m = ⌊ (a n + b) / c ⌋$ 个 $U$ 组成；
第 $i$ 个 $R$ 前方的 $U$ 的数量恰等于 $⌊ (a i + b) / c ⌋$ ，其中， $i = 1, \dots, n$ 。

从几何直观上看，这大致相当于从原点开始，每向右穿过一次竖向的网格线，就写下一个 $R$ ，每向上穿过一次横向的网格线，就写下一个 $U$ 。如下图所示：

当然，这样的定义还需要考量一系列特殊情形：

经过整点（即同时上穿和右穿）时，需要先写 $U$ 再写 $R$ ；
字符串开始时，除了在 $(0, 1]$ 区间内上穿网格线的次数外，还需要格外补充 $⌊ b / c ⌋$ 个 $U$ ；
字符串结束时，不能有格外的 $U$ 。

如果对于几何直观的描述有任何不明晰的地方，可以参考上述代数方法的定义辅助理解。几何直观的描述，有助于理解下文的算法过程。

万能欧几里得算法的基本思路，就是将操作序列中的 $U$ 和 $R$ 都视作某个幺半群内的元素，将整个操作序列视为幺半群内元素的乘积，而问题最终的答案与这个乘积有关。

比如，本题中，可以定义状态向量 $v = (1, y, \sum y)$ ，表示自原点开始，经历了若干次上穿和右穿网格线后，当前的状态。其中，第一个分量是常数，第二个分量是纵坐标 $y$ ，第三个分量是要求的和式。起始时，有 $v = (1, 0, 0)$ 。每向上穿过一次网格线，纵坐标就累加一，即相当于将状态向量右乘以矩阵

U = (\begin{matrix} 1 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) .

每向右穿过一次网格线，和式就累加一次纵坐标，即相当于将状态向量右乘以矩阵

R = (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{matrix}) .

因此，最终的状态就是乘积 $(1, 0, 0) S$ ，其中， $S$ 理解为上述矩阵的乘积。所求的答案，就是最终状态的第三个分量。

除了将幺半群中的元素定义为矩阵以外，还可以将它们定义为一段操作序列对于最终结果的贡献，然后将操作的乘积定义为两段操作序列的贡献的合并。

本题中，可以定义每段操作序列的贡献为 $(x, y, \sum y)$ 。为了严谨地解释这些记号，可以将这些分量都看作是操作序列的函数，也就是说，对于操作序列 $S$ ，它的贡献可以写作 $(x (S), y (S), (\sum y) (S))$ 。其中， $x (S)$ 和 $y (S)$ 分别对应着操作序列 $S$ 中 $R$ 和 $U$ 的数量，也就是该线段右穿和上穿网格线的次数。最后一项中的求和符号，一般地，有如下定义：对于操作序列上的函数 $f (S)$ ，可以定义 $(\sum f) (S)$ ，或记作 $\sum_{S} f$ ，为下面的表达式：

\sum_{S} f := \sum {f (S_{[1, r]}) : S_{r} = R} .

其中， $S_{r}$ 是 $S$ 中的第 $r$ 个字符， $S_{[1, r]}$ 是 $S$ 中前 $r$ 个字符组成的前缀。也就是说，这个求和记号，可以看作是对于操作序列 $S$ 中所有以 $R$ 结尾的前缀进行的求和。比如说，有

\sum_{S} 1 = x, \sum_{S} x = \frac{1}{2} x (x + 1) .

再比如说， $\sum y$ 就是操作序列中，每次右穿网格线时，之前上穿网格线的次数的累加。对于整段操作序列来说， $y$ 在所有以 $R$ 结尾的前缀处的值，就是在 $i = 1, \dots, n$ 处的所有 $⌊ (a i + b) / c ⌋$ 的值。因此，对于整段操作序列计算的 $\sum y$ ，就是本题最终要求的量。

初始时，有 $U = (0, 1, 0)$ ， $R = (1, 0, 0)$ 。进一步，可以将两个元素 $(x_{1}, y_{1}, s_{1})$ 和 $(x_{2}, y_{2}, s_{2})$ 的乘积定义为

(x_{1}, y_{1}, s_{1}) \cdot (x_{2}, y_{2}, s_{2}) = (x_{1} + x_{2}, y_{1} + y_{2}, s_{1} + s_{2} + x_{2} y_{1}) .

其中，最后一项贡献合并的结果可以通过如下计算得到：

\sum_{S_{1} + S_{2}} y = \sum_{S_{1}} y + \sum_{S_{2}} (y + y_{1}) = \sum_{S_{1}} y + \sum_{S_{2}} y + y_{1} \sum_{S_{2}} 1 = s_{1} + s_{2} + x_{2} y_{1} .

容易验证，这个乘法运算满足结合律，且幺元为 $(0, 0, 0)$ ，所以这些元素在该乘法运算下构成幺半群。所求的答案，就是乘积的第三个分量。

这两种方法都可以得到正确的结果。但是，因为保留了较多的冗余信息，矩阵运算的常数较大，所以第二种方法在处理实际问题时更为实用。

算法过程

与类欧几里得算法整体约化不同，万能欧几里得算法约化问题的手段是将这些操作分批次地合并。记字符串对应的操作的乘积为

F (a, b, c, n, U, R) .

约化过程具体如下：

当 $b \geq c$ 时，操作序列的开始有 $⌊ b / c ⌋$ 个 $U$ ，直接计算它们的乘积，并将这些 $U$ 从操作序列中移除。此时，第 $i$ 个 $R$ 前方的 $U$ 的数量等于
$⌊ \frac{a i + b}{c} ⌋ - ⌊ \frac{b}{c} ⌋ = ⌊ \frac{a i + (b mod c)}{c} ⌋ .$
因此，这相当于将线段参数由 $(a, b, c, n)$ 变为 $(a, b mod c, c, n)$ 。所以，对于这种情形，有
$F (a, b, c, n, U, R) = U^{⌊ b / c ⌋} F (a, b mod c, c, n, U, R) .$
当 $a \geq c$ 时，操作序列中每个 $R$ 的前方都至少有 $⌊ a / c ⌋$ 个 $U$ ，可以将它们合并到 $R$ 上。也就是说，可以用 $U^{⌊ a / c ⌋} R$ 替代 $R$ 。合并后的字符串中，第 $i$ 个 $R$ 前方的 $U$ 的数量等于
$⌊ \frac{a i + b}{c} ⌋ - ⌊ \frac{a}{c} ⌋ i = ⌊ \frac{(a mod c) i + b}{c} ⌋ .$
因此，这相当于将线段参数由 $(a, b, c, n)$ 变为 $(a mod c, b, c, n)$ 。所以，对于这种情形，有
$F (a, b, c, n, U, R) = F (a mod c, b, c, n, U, U^{⌊ a / c ⌋} R) .$
对于剩下的情形，需要翻转横纵坐标，这基本是在交换 $U$ 和 $R$ ，只是翻转后线段的参数需要仔细计算。结合操作序列的定义可知，需要确定系数 $(a^{'}, b^{'}, c^{'}, n^{'})$ 使得变换前的操作序列中，第 $j$ 个 $U$ 前方的 $R$ 的数量恰为 $⌊ (a^{'} j + b^{'}) / c^{'} ⌋$ 且总共有 $n^{'}$ 个 $U$ 。根据定义可知，
$n^{'} = ⌊ \frac{a n + b}{c} ⌋ = m,$ $n'=\left\lfloor\dfrac{an+b}{c}\right\rfloor = m,$
而第 $j$ 个 $U$ 前方的 $R$ 的数量，就等于最大的 $i$ 使得
$\begin{aligned} ⌊ \frac{a i + b}{c} ⌋ < j & ⟺ \frac{a i + b}{c} < j ⟺ i < \frac{c j - b}{a} \\ ⟺ i < ⌈ \frac{c j - b}{a} ⌉ = ⌊ \frac{c j - b - 1}{a} ⌋ + 1. \end{aligned}$ $\begin{aligned} \left\lfloor\dfrac{ai+b}{c}\right\rfloor < j &\iff \dfrac{ai+b}{c} < j \iff i < \dfrac{cj-b}{a} \\ &\iff i < \left\lceil\dfrac{cj-b}{a}\right\rceil = \left\lfloor\dfrac{cj-b - 1}{a}\right\rfloor + 1. \end{aligned}$
因此， $i = ⌊ (c j - b - 1) / a ⌋$ 。这一推导过程与前文类欧几里得算法的推导类似，同样利用了上下取整函数的性质。

有两处细节需要处理：
- 截距项 $- (b + 1) / a$ 为负数。注意到，如果将线段向左平移一个单位，就可以让截距项恢复为非负数，因为总有 $(c - b - 1) / a \geq 0$ 。因此，可以将交换前的第一段 $R^{⌊ (c - b - 1) / a ⌋} U$ 提取出来，只交换剩余操作序列中的 $U$ 和 $R$ ；
- 交换 $U$ 和 $R$ 后，结尾存在多余的 $U$ 。因此，交换 $U$ 和 $R$ 之前，需要首先将最后一段 $R$ 提取出来，只交换剩余操作序列中的 $U$ 和 $R$ 。这一段 $R$ 的数量为 $n - ⌊ (c m - b - 1) / a ⌋$ 。
去掉头尾若干个字符后，第 $j$ 个 $U$ 前方的 $R$ 的数量变为：
$⌊ \frac{c (j + 1) - b - 1}{a} ⌋ - ⌊ \frac{c - b - 1}{a} ⌋ = ⌊ \frac{c j + (c - b - 1) mod a}{a} ⌋ .$ $\left\lfloor\dfrac{c(j+1)-b-1}{a}\right\rfloor - \left\lfloor\dfrac{c-b-1}{a}\right\rfloor = \left\lfloor\dfrac{cj+(c-b-1)\bmod a}{a}\right\rfloor.$
回忆起，交换前的序列中 $U$ 的数量为 $m = ⌊ (a n + b) / c ⌋$ 。而上述左移一个单位的操作，要求保证交换前至少存在一个 $U$ ，也就是 $m > 0$ 。利用这一条件，可以分为两种情形：
- 对于 $m > 0$ 的情形，处理了上面的两点后，交换完 $U$ 和 $R$ 的操作序列就是对应着参数为 $(c, (c - b - 1) mod a, a, m - 1)$ 的线段的合法序列。所以，有
  $F (a, b, c, n, U, R) = R^{⌊ (c - b - 1) / a ⌋} U F (c, (c - b - 1) mod a, a, m - 1, R, U) R^{n - ⌊ (c m - b - 1) / a ⌋} .$
- 特别地，对于 $m = 0$ 的情形，交换前的操作序列中只包含 $n$ 个 $R$ ，无需交换，可以直接返回：
  $F (a, b, c, n, U, R) = R^{n} .$
  与类欧几里得算法不同，万能欧几里得算法的这一特殊情形需要单独处理，否则会因涉及负幂次而无法正确计算。

利用这些讨论，就可以将问题递归地解决。

假设幺半群内元素单次相乘的时间复杂度是 $O (1)$ 的。那么，如果计算过程中这些元素的幂次计算都使用快速幂进行，最终的算法复杂度就是 $O (\log max {a, c} + \log (b / c))$ 的[^complexity]。

对复杂度的解释

对比（类）欧几里得算法，万能欧几里得算法只是多了求快速幂的步骤。其余的计算过程的复杂度和类欧几里得算法相仿，已经说明是 $O (\log min {a, c, n})$ 的。现在，需要计算这些快速幂的总复杂度。

除了第一轮迭代，都有 $b

本页面最近更新：2025/10/13 16:54:57，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：sshwy, StudyingFather, Enter-tainer, H-J-Granger, Tiphereth-A, countercurrent-time, NachtgeistW, Early0v0, Ir1d, MegaOwIer, Xeonacid, AngelKitty, c-forrest, CCXXXI, cjsoft, diauweb, ezoixx130, FFjet, GekkaSaori, Henry-ZHR, Konano, LovelyBuggies, Makkiy, mgt, minghu6, P-Y-Y, PotassiumWings, SamZhangQingChuan, Suyun514, weiyong1024, alphagocc, cxm1024, GavinZhengOI, Gesrua, Great-designer, iamtwz, ksyx, kxccc, lychees, megakite, Peanut-Tang, r-value, SukkaW, TonyYin0418
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用