类欧几里德算法

引入

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

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

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

类欧几里得算法

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

其中， 都是正整数。

代数解法

首先，将 对 取模，可以简化问题，将问题转化为 的情形：

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

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

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

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

令 ，这就又回到了前面讨论过的 的情形。

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

在计算过程中，可能会出现 的情形，此时内层递归会出现 。这并不影响最终的结果。但是，如果要求出现 时，直接终止算法，那么算法的时间复杂度可以改良为 的。

对复杂度的解释

利用该算法和欧几里得算法的相似性，很容易说明它的时间复杂度是 的。因此，只需要说明，如果在 时终止算法，那么它的时间复杂度也是 的。

令 ，并记 ，，它们分别相当于几何直观（见下一节）中点阵图的面积和直线的斜率。对于充分大的 ，近似有 。

考察 和 在算法过程中的变化。第一步取模时， 保持不变， 近似由 变为 ，相当于斜率由 变为 ，而 也近似变为原来的 倍。第二步交换横纵坐标时， 近似保持不变， 则变为它的倒数。因此，若设两步操作后，二元组 变为 ，则有 且 。

因为 ，所以，递归计算两轮后，乘积缩小的倍数最少为

因此，至多 轮，算法必然终止。因为从第二轮开始，每轮开始时的 总是不超过上一轮取模结束后的 ，而后者大致为 且 ，因而 。这就得到了上述结论。

模板题的参考实现如下：

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

#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;
}

几何直观

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

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

下方， 轴上方（不包括 轴），且横坐标位于 之间的格点数目。

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

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

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

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

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

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

这也是斜率为有理数时的类欧几里得算法的复杂度是 的原因。

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

例题

【模板】类欧几里得算法

多组询问。给定正整数 ，求

解答一

类似于 的推导，可以得到 的递归表达式。

首先，利用取模，将问题转化为 的情形：

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

那么，对于和式 ，有

对于和式 ，有

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

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

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

#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

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

解答一

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

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

其中的函数 具有形式

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

同样分为两种情形讨论。如果 ，那么

问题转化为斜率小于一的情形。如果 ，那么设 ，有

此处的推导中，交换 和 的条件比正文中的情形更为简单，是因为直线 上没有除了原点之外的格点。关键在于交换后的求和式写成 的形式，这相当于要求 满足

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

因此，有

这说明

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

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

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

#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; }