企业官网建设流程全解析

类欧几里得算法与万能欧几里得算法：OI-wiki 中直线下整点计数与操作序列方法

【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. （某大型游戏线上攻略，内含炫酷算术魔法）项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki

本文围绕 OI-wiki 数论模块中的 euclidean.md 展开，系统讲解用于计算 $\left\lfloor\dfrac{ai+b}{c}\right\rfloor$ 形式求和的类欧几里得算法，以及将其推广、可求解更多带权和式的万能欧几里得算法。这两类算法利用分数自身的递归结构与欧几里得算法的内在联系，将大范围求和问题约化为 $O(\log\min{a,c,n})$ 的递归过程，是竞赛编程中处理「直线下方整点计数」类问题的利器。读完本文，你将掌握 $f/g/h$ 三类和式的代数推导、几何直观理解，以及基于幺半群抽象出的统一递归模板，并能在 Library Checker - Sum of Floor of Linear 与 Luogu P5170 等模板题上直接落地实现。

引入：从欧几里得算法到类欧几里得算法

类欧几里得算法由洪华敦在 2016 年冬令营营员交流中提出，常用于解决形如

$$ \left\lfloor\dfrac{ai+b}{c}\right\rfloor $$

结构的数列（下标为 $i$）的求和问题。它的核心思想是：利用分数自身的递归结构，将问题转化为更小规模的问题递归求解。

之所以冠以「类欧几里得」之名，是因为分数的递归结构与 欧几里得算法 存在直接联系（详见 连分数表示的求法）。事实上，连分数 和 Stern–Brocot 树 等方法同样刻画了分数的递归结构，因此能用类欧几里得算法解决的问题通常也可以用这些方法解决；但相较之下，类欧几里得算法通常更容易理解，实现也更为简明。

类欧几里得算法

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

$$ f(a,b,c,n)=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor, $$

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

代数解法

第一步：取模约化。将 $a,b$ 对 $c$ 取模，可以简化问题，将问题转化为 $0\le a,b<c$ 的情形：

$$ \begin{aligned} f(a,b,c,n)&=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor\ &=\sum_{i=0}^n\left\lfloor \frac{\left(\left\lfloor\frac{a}{c}\right\rfloor c+(a\bmod c)\right)i+\left(\left\lfloor\frac{b}{c}\right\rfloor c+(b\bmod c)\right)}{c}\right\rfloor\ &=\sum_{i=0}^n\left(\left\lfloor\frac{a}{c}\right\rfloor i+\left\lfloor\frac{b}{c}\right\rfloor+\left\lfloor\frac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c} \right\rfloor\right)\ &=\frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor +(n+1)\left\lfloor\frac{b}{c}\right\rfloor+f(a\bmod c,b\bmod c,c,n). \end{aligned} $$

其中 $\sum_{i=0}^n i=\frac{n(n+1)}{2}$ 的等差数列求和是唯一的「闭合解」部分。

第二步：交换求和次序。现在考虑转化后 $0\le a,b<c$ 的问题。令

$$ m = \left\lfloor \frac{an+b}{c} \right\rfloor. $$

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

$$ \sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor =\sum_{i=0}^n\sum_{j=0}^{m-1}\left[j<\left\lfloor \frac{ai+b}{c} \right\rfloor\right]. $$

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

$$ \begin{aligned} &j<\left\lfloor \frac{ai+b}{c} \right\rfloor = \left\lceil \frac{ai+b+1}{c} \right\rceil-1\ &\iff j + 1 < \left\lceil \frac{ai+b+1}{c} \right\rceil \iff j+1< \frac{ai+b+1}{c} \ &\iff \dfrac{cj+c-b-1}{a} < i \iff \left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor < i. \end{aligned} $$

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

$$ \begin{aligned} f(a,b,c,n)&=\sum_{j=0}^{m-1} \sum_{i=0}^n\left[i>\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor \right]\ &=\sum_{j=0}^{m-1}\left(n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)\ &=nm-f\left(c,c-b-1,a,m-1\right). \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$ 的边界情形。计算过程中可能出现 $m=0$，此时内层递归会出现 $n=-1$，但这不影响最终结果。如果要求出现 $m=0$ 时直接终止算法，算法的时间复杂度可以改良为 $O(\log\min{a,c,n})$。

复杂度的几何解释。利用该算法与欧几里得算法的相似性，容易说明其时间复杂度是 $O(\log\min{a,c})$；而若在 $m=0$ 时终止算法，还需说明它也是 $O(\log n)$ 的。令 $m=\lfloor(an+b)/c\rfloor$，记 $S=mn$、$k=m/n$，它们分别相当于几何直观中（见下一小节）点阵图的面积和直线的斜率；对于充分大的 $n$，近似有 $k\doteq a/c$。

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

因为 $1\le\lfloor k'\rfloor\le k'<\lfloor k'\rfloor+1$，递归计算两轮后乘积缩小的倍数最少为

$$ (k'-\lfloor k'\rfloor)(k-\lfloor k'\rfloor) = 1-\dfrac{\lfloor k'\rfloor}{k'} < 1-\dfrac{\lfloor k'\rfloor}{\lfloor k'\rfloor+1} = \dfrac{1}{\lfloor k'\rfloor+1}\le \dfrac{1}{2}. $$

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

模板题参考实现。仓库中的 euclidean-0.cpp 是在 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; }

注意主函数中读取顺序为n c a b且传入n - 1，这是 Library Checker 题目中求和范围 $[0,n)$ 的约定（原题下标从 $0$ 到 $n-1$）。

几何直观

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

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

$$ y = \dfrac{ax+b}{c} $$

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

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

$$ f(a,b,c,n) = f(a\bmod c,b\bmod c,c,n) + \dfrac{1}{2}(n+1)\left(\left\lfloor\dfrac{b}{c}\right\rfloor+\left(\left\lfloor\dfrac{a}{c}\right\rfloor n+\left\lfloor\dfrac{b}{c}\right\rfloor\right)\right). $$

第二步：翻转横纵坐标轴。如下图最左部分所示，红点和蓝点构成一个横向长度为 $n$、纵向长度为 $m=\lfloor(an+b)/c\rfloor$ 的矩形点阵。要计算红点数量，只需计算蓝点数量再用矩形点阵总数减去即可。翻转后，左半部分的蓝点点阵变成某条直线下方的红色点阵；且翻转后斜率大于一，又回到上文已处理的情形。

关键在于新红色点阵上方直线的方程。将最左部分的横纵坐标轴翻转得到中间部分：翻转后的红色点阵上方的直线（中间部分实线）并非翻转前直线（最左部分实线）的直接翻转，而是向左上平移一点点的结果（最左部分虚线）。这是因为直接将直线翻转会得到中间部分虚线，而按定义它下方的格点包含恰好落在直线上的格点，会造成重复计数。为避免这一点，需将翻转后得到的直线 $y=(cx-b)/a$ 向下平移一点点得到 $y=(cx-b-1)/a$，这样它下方的点阵才恰为翻转前的蓝色点阵。

还有一处细节：中间部分直线的截距是负数，尚未回到初始情形。要让截距恢复非负，只需将直线向左平移一个单位——这不会漏掉任何格点，因为翻转前的蓝色点阵中没有纵坐标为零的点，翻转后也就不存在横坐标为零的点。最终直线方程变为 $y=(cx+c-b-1)/a$，点阵横坐标上界也从 $m$ 变为 $m-1$。这一步骤归纳为：

$$ f(a,b,c,n) = mn - 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})$ 的原因。利用类似的几何直观，还可以将类欧几里得算法推广到斜率为无理数的情形（见后文例题）。

例题一：【模板】类欧几里得算法（Luogu P5170）

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

$$ \begin{aligned} f(a,b,c,n) &= \sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor,\ g(a,b,c,n) &= \sum_{i=0}^ni\left\lfloor \frac{ai+b}{c} \right\rfloor,\ h(a,b,c,n) &= \sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor^2. \end{aligned} $$

推导 $g,h$ 的递归表达式。类似于 $f$ 的推导，首先利用取模将问题转化为 $0\le a,b<c$ 的情形：

$$ \begin{aligned} g(a,b,c,n) &=g(a\bmod c,b\bmod c,c,n)+\left\lfloor\frac{a}{c}\right\rfloor\frac{n(n+1)(2n+1)}{6}+\left\lfloor\frac{b}{c}\right\rfloor\frac{n(n+1)}{2}, \ h(a,b,c,n)&=h(a\bmod c,b\bmod c,c,n)\ &\quad+2\left\lfloor\frac{b}{c}\right\rfloor f(a\bmod c,b\bmod c,c,n) +2\left\lfloor\frac{a}{c}\right\rfloor g(a\bmod c,b\bmod c,c,n)\ &\quad+\left\lfloor\frac{a}{c}\right\rfloor^2\frac{n(n+1)(2n+1)}{6}+\left\lfloor\frac{b}{c}\right\rfloor^2(n+1) +\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor n(n+1). \end{aligned} $$

然后利用交换求和次序进一步转化。同样令 $m = \left\lfloor \frac{an+b}{c} \right\rfloor$。对于 $g$：

$$ \begin{aligned} g(a,b,c,n)&=\sum_{i=0}^ni\left\lfloor \frac{ai+b}{c} \right\rfloor\ &=\sum_{i=0}^n \sum_{j=0}^{m-1}i \left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor\right] \ &=\sum_{j=0}^{m-1}\sum_{i=0}^n i\left[i>\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor \right]\ &=\sum_{j=0}^{m-1}\dfrac{1}{2}\left(\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor+n+1\right)\left(n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)\ &=\dfrac{1}{2}mn(n+1) - \dfrac{1}{2}\sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor - \dfrac{1}{2}\sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor^2\ &=\dfrac{1}{2}mn(n+1) - \dfrac{1}{2}f(c,c-b-1,a,m-1) - \dfrac{1}{2}h(c,c-b-1,a,m-1). \end{aligned} $$

对于 $h$：

$$ \begin{aligned} h(a,b,c,n)&=\sum_{i=0}^n\left\lfloor \frac{ai+b}{c} \right\rfloor^2\ &=\sum_{i=0}^n\sum_{j=0}^{m-1}(2j+1)\left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor\right]\ &=\sum_{j=0}^{m-1}\sum_{i=0}^n(2j+1)\left[i>\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor \right]\ &=\sum_{j=0}^{m-1}(2j+1)\left(n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)\ &=nm^2 - \sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor - 2\sum_{j=0}^{m-1}j\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\ &=nm^2 - f(c,c-b-1,a,m-1) - 2g(c,c-b-1,a,m-1). \end{aligned} $$

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

$$ \sum_{i=0}^ni^r\left\lfloor \frac{ai+b}{c} \right\rfloor^s = \sum_{i=0}^n\sum_{j=0}^{m-1} i^r\left((j+1)^s-j^s\right)\left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor\right]. $$

本题的另一个特点是 $g$ 和 $h$ 在递归计算时相互交错，因此需要将 $(f,g,h)$ 作为三元组同时递归。仓库中的 euclidean-1.cpp 实现了这一做法，它在模 $998244353$ 意义下计算（使用 $i2=(M+1)/2$、$i6=(M+1)/6$ 处理逆元），并在取模约化分支中显式叠加了 $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; }

注意 $g$ 分支中m * n2 - (tmp.f + tmp.h)/2与推导式对应，$h$ 分支中n*m*m - tmp.f - 2*tmp.g与推导式对应，输出顺序为f h g（题目要求）。

例题二：【清华集训 2014】Sum（Luogu P5172）

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

$$ \sum_{d=1}^n(-1)^{\lfloor d\sqrt{r}\rfloor}. $$

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

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

$$ \begin{aligned} \sum_{d=1}^n(-1)^{\lfloor d\sqrt{r}\rfloor} &= \sum_{d=1}^n\left(1 - 2(\lfloor d\sqrt{r}\rfloor\bmod 2)\right)\ &= n - 2\sum_{d=1}^n\left(\lfloor d\sqrt{r}\rfloor - 2\left\lfloor\dfrac{\lfloor d\sqrt{r}\rfloor}{2}\right\rfloor\right) \ &= n - 2\sum_{d=1}^n\lfloor d\sqrt{r}\rfloor + 4 \sum_{d=1}^n\left\lfloor\dfrac{d\sqrt{r}}{2}\right\rfloor\ &= n - 2f(n,1,0,1) + 4f(n,1,0,2) \end{aligned} $$

其中函数 $f$ 具有形式

$$ f(a,b,c,n) = \sum_{i=1}^n\left\lfloor\dfrac{a\sqrt{r}+b}{c}i\right\rfloor. $$

与正文算法不同，此处斜率不再是有理数。设斜率 $k = \dfrac{a\sqrt{r}+b}{c}$，分两种情形讨论。

若 $k\ge 1$：

$$ \begin{aligned} f(a,b,c,n) &= \sum_{i=1}^n \lfloor ki\rfloor = \sum_{i=1}^n \lfloor(k-\lfloor k\rfloor)i\rfloor + \lfloor k\rfloor \sum_{i=1}^ni\ &= \lfloor k\rfloor\dfrac{n(n+1)}{2} + f(a,b-c\lfloor k\rfloor,c,n). \end{aligned} $$

问题转化为斜率小于一的情形。

若 $k<1$：设 $m=\lfloor nk\rfloor$，有

$$ \begin{aligned} f(a,b,c,n) &= \sum_{i=1}^n \lfloor ki\rfloor = \sum_{i=1}^n\sum_{j=1}^m[j\le\lfloor ki\rfloor]\ &= \sum_{j=1}^m\sum_{i=1}^n[i>\lfloor k^{-1}j\rfloor] = nm - \sum_{j=1}^m\sum_{i=1}^n[i\le\lfloor k^{-1}j\rfloor]. \end{aligned} $$

此处交换 $i,j$ 的条件比正文更简单，是因为直线 $y=kx$ 上除原点外没有格点。关键在于将交换后的求和式写成 $f(a,b,c,n)$ 的形式，即要求 $a',b',c'$ 满足

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

分母有理化即可得到

$$ k^{-1} = \dfrac{c}{a\sqrt{r}+b} = \dfrac{ca\sqrt{r}-cb}{a^2r-b^2}. $$

因此有

$$ a'=ca,~b'=-cb,~c'=a^2r-b^2, $$

即

$$ f(a,b,c,n) = nm - f(ca,-cb,a^2r-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)$ 的。仓库中的 euclidean-2.cpp 实现了上述过程，其solve中对完全平方数单独分支、f中先约分 $\gcd(a,b,c)$、用sqrtl计算斜率并据此分支取整，读者可对照阅读。

例题三：Fraction（Luogu P5179）

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

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

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

$$ \left\lfloor\dfrac{a}{b}\right\rfloor \le \dfrac{a}{b} <\dfrac{p}{q} <\dfrac{c}{d}\le\left\lfloor\dfrac{a}{b}\right\rfloor+1. $$

从这个不等式可以看出 $p/q$ 的整数部分可确定为 $\lfloor a/b\rfloor$，直接消去整数部分后整体取倒数，用于确定小数部分——这正是确定 $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$。但这说明

$$ \dfrac{a}{b} < \dfrac{r}{s} < \dfrac{r}{q} < \dfrac{p}{q} < \dfrac{c}{d}, $$

即 $r/q$ 无论按哪个字典序都严格更小，与所设矛盾。因此上述算法是正确的。仓库中的 euclidean-3.cpp 给出了极简实现，递归中交换 $p,q$ 与 $a,b$ 的角色并在回溯时恢复整数部分。

万能欧几里得算法

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

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

$$ f(a,b,c,n)=\sum_{i=1}^n\left\lfloor \frac{ai+b}{c} \right\rfloor, $$

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

问题转化：操作序列与幺半群

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

$$ y = \frac{ax+b}{c},~0< x\le n. $$

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

• 字符串恰有 $n$ 个 $R$ 和 $m=\lfloor(an+b)/c\rfloor$ 个 $U$；
• 第 $i$ 个 $R$ 前方的 $U$ 的数量恰等于 $\lfloor(ai+b)/c\rfloor$，其中 $i=1,\cdots,n$。

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

这样定义还需考量一系列特殊情形：

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

如果几何直观描述有不明晰之处，可参考上述代数方法的定义辅助理解。

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

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

$$ U = \begin{pmatrix}1 & 1 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1\end{pmatrix}. $$

每向右穿过一次，和式累加一次纵坐标，相当于右乘矩阵

$$ R = \begin{pmatrix}1 & 0 & 0 \ 0 & 1 & 1 \ 0 & 0 & 1\end{pmatrix}. $$

最终状态即乘积 $(1,0,0)S$（$S$ 为上述矩阵的乘积），所求答案就是最终状态的第三个分量。

方式二：贡献合并（更实用）。除了矩阵，还可以将幺半群元素定义为一段操作序列对最终结果的贡献，将操作乘积定义为两段贡献的合并。本题中定义每段操作序列的贡献为 $(x,y,\sum y)$，其中 $x(S)$、$y(S)$ 分别对应 $S$ 中 $R$ 和 $U$ 的数量，最后一项的求和符号一般定义如下：对于操作序列上的函数 $f(S)$，定义

$$ \sum_S f := \sum{f(S_{[1,r]}):S_r=R}, $$

其中 $S_r$ 是 $S$ 中第 $r$ 个字符，$S_{[1,r]}$ 是前 $r$ 个字符组成的前缀，即对操作序列中所有以 $R$ 结尾的前缀求和。例如

$$ \sum_S 1 = x,~ \sum_S x = \dfrac{1}{2}x(x+1). $$

而 $\sum y$ 就是每次右穿时之前上穿次数的累加：对于整段操作序列，$y$ 在所有以 $R$ 结尾的前缀处的值，正是 $i=1,\cdots,n$ 处的所有 $\lfloor(ai+b)/c\rfloor$ 值，因此整段序列的 $\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_2y_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_2y_1 $$

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

两种方法都能得到正确结果，但矩阵运算保留了较多冗余信息、常数较大，因此第二种方法（贡献合并）在实际问题中更为实用。

算法过程：分批次合并操作

与类欧几里得算法整体约化不同，万能欧几里得算法约化问题的手段是将操作分批次合并。记字符串对应的操作的乘积为 $F(a,b,c,n,U,R)$，约化过程如下：

情形一：$b\ge c$。操作序列开始有 $\lfloor b/c\rfloor$ 个 $U$，直接计算其乘积并移除。此时第 $i$ 个 $R$ 前方的 $U$ 数量等于

$$ \left\lfloor\dfrac{ai+b}{c}\right\rfloor - \left\lfloor\dfrac{b}{c}\right\rfloor = \left\lfloor\dfrac{ai+(b\bmod c)}{c}\right\rfloor, $$

相当于线段参数由 $(a,b,c,n)$ 变为 $(a,b\bmod c,c,n)$。因此

$$ F(a,b,c,n,U,R) = U^{\lfloor b/c\rfloor}F(a,b\bmod c,c,n,U,R). $$

情形二：$a\ge c$。每个 $R$ 前方都至少有 $\lfloor a/c\rfloor$ 个 $U$，可将其合并到 $R$ 上，即用 $U^{\lfloor a/c\rfloor}R$ 替代 $R$。合并后第 $i$ 个 $R$ 前方的 $U$ 数量等于

$$ \left\lfloor\dfrac{ai+b}{c}\right\rfloor - \left\lfloor\dfrac{a}{c}\right\rfloor i = \left\lfloor\dfrac{(a\bmod c)i+b}{c}\right\rfloor, $$

相当于参数由 $(a,b,c,n)$ 变为 $(a\bmod c,b,c,n)$。因此

$$ F(a,b,c,n,U,R) = F(a\bmod c,b,c,n,U,U^{\lfloor a/c\rfloor}R). $$

情形三：其余情形（翻转横纵坐标）。这基本是在交换 $U$ 和 $R$，但翻转后的参数需要仔细计算。结合操作序列定义，需确定系数 $(a',b',c',n')$ 使变换前的操作序列中第 $j$ 个 $U$ 前方的 $R$ 数量恰为 $\lfloor(a'j+b')/c'\rfloor$ 且总共有 $n'$ 个 $U$。根据定义

$$ n'=\left\lfloor\dfrac{an+b}{c}\right\rfloor = m, $$

而第 $j$ 个 $U$ 前方的 $R$ 数量等于最大的 $i$ 使得

$$ \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 = \lfloor(cj-b-1)/a\rfloor$。这一推导与前文类欧几里得算法类似，同样利用了上下取整函数的性质。

有两处细节需要处理：

• 负截距。截距项 $-(b+1)/a$ 为负数。注意到将线段向左平移一个单位可使截距恢复非负，因为总有 $(c-b-1)/a\ge 0$。因此可将交换前的第一段 $R^{\lfloor(c-b-1)/a\rfloor}U$ 提取出来，只交换剩余操作序列中的 $U$ 和 $R$；
• 结尾多余的 $U$。交换 $U$ 和 $R$ 后结尾存在多余的 $U$，因此交换前需先将最后一段 $R$ 提取出来（其数量为 $n-\lfloor(cm-b-1)/a\rfloor$），只交换剩余部分。

去掉头尾若干字符后，第 $j$ 个 $U$ 前方的 $R$ 数量变为：

$$ \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=\lfloor(an+b)/c\rfloor$，而左移操作要求交换前至少存在一个 $U$，即 $m>0$。据此分两种情形：

• $m>0$：处理上述两点后，交换完 $U$ 和 $R$ 的操作序列就是参数为 $(c,(c-b-1)\bmod a,a,m-1)$ 的线段的合法序列，所以

$$ F(a,b,c,n,U,R) = R^{\lfloor(c-b-1)/a\rfloor}UF(c,(c-b-1)\bmod a,a,m-1,R,U)R^{n-\lfloor(cm-b-1)/a\rfloor}. $$

• $m=0$：交换前的操作序列只包含 $n$ 个 $R$，无需交换，直接返回

$$ F(a,b,c,n,U,R) = R^n. $$

与类欧几里得算法不同，万能欧几里得算法的这一特殊情形需要单独处理，否则会因涉及负幂次而无法正确计算。

利用这些讨论即可递归求解。

复杂度。假设幺半群内元素单次相乘为 $O(1)$，且元素幂次计算都使用 快速幂，最终算法复杂度为 $O(\log\max{a,c}+\log(b/c))$。

复杂度论证的要点是：除第一轮迭代外都有 $b<c$，每轮迭代涉及三次快速幂，其总复杂度为

$$ O\left(\log\left\lfloor\dfrac{a}{c}\right\rfloor+\log\left\lfloor\dfrac{c-b_1-1}{a_1}\right\rfloor+\log\left(n-\left\lfloor\dfrac{cm-b_1-1}{a_1}\right\rfloor\right)\right), $$

其中 $a_1=a\bmod c$、$b_1=b\bmod c$ 且 $m=\lfloor(a_1n+b_1)/c\rfloor$。后两项分别有估计

$$ \begin{aligned} \dfrac{c-b_1-1}{a_1} &\le \dfrac{c}{a_1},\ n-\left\lfloor\dfrac{cm-b_1-1}{a_1}\right\rfloor &\le n - \dfrac{cm-b_1-1}{a_1} + 1 \ &\le n - \dfrac{c((a_1n+b_1)/c-1)-b_1-1}{a_1} +1 \ &= \dfrac{c+1}{a_1}+1, \end{aligned} $$

故这两项复杂度都是 $O(\log(c/a_1))$。每一轮迭代中线段参数由 $(a,\cdot,c,\cdot)$ 变换为 $(c,\cdot,a\bmod c,\cdot)$，该轮总时间复杂度为

$$ O\left(\log\dfrac{a}{c}+\log\dfrac{c}{a\bmod c}\right), $$

全部递归轮次中这些项可以裂项相消，总和为 $O(\log a+\log c)=O(\log\max{a,c})$。再加上第一轮迭代中 $U^{\lfloor b/c\rfloor}$ 的快速幂复杂度 $O(\log(b/c))$，即得总复杂度 $O(\log\max{a,c}+\log(b/c))$。

关于 $O(\log(b/c))$ 项的注记。通常考虑的问题中 $b$ 与 $a$ 同阶，这一项可以忽略；而且如果在调用万能欧几里得算法前先进行一轮类欧几里得算法的取模消除 $b$ 的影响，该项快速幂的复杂度可以规避。这其实是因为通常问题中 $U$ 的初始形式较为特殊，其幂次有更简单的形式，不需要通过快速幂计算——比如正文例子中 $U^{\lfloor b/a\rfloor}$ 的结果，就是将 $U$ 中不在对角线上的那个 $1$ 替换为 $\lfloor b/a\rfloor$，无需快速幂。

统一模板。万能欧几里得算法的流程可以写成统一模板，处理具体问题时只需更改模板类型T的实现。仓库中的 euclidean-4.cpp 完整实现了这一模板，其核心euclid函数（文件内标记为euclidean片段）与上述三种情形一一对应：

// Class T implements the monoid. // Assume that it provides a multiplication operator // and a default constructor returning the unity in the monoid. // Binary exponentiation. template <typename T> T pow(T a, int b) { T res; for (; b; b >>= 1) { if (b & 1) res = res * a; a = a * a; } return res; } // Universal Euclidean algorithm. template <typename T> T euclid(int a, int b, int c, int n, T U, T R) { if (b >= c) return pow(U, b / c) * euclid(a, b % c, c, n, U, R); if (a >= c) return euclid(a % c, b, c, n, U, pow(U, a / c) * R); auto m = ((long long)a * n + b) / c; if (!m) return pow(R, n); return pow(R, (c - b - 1) / a) * U * euclid(c, (c - b - 1) % a, a, m - 1, R, U) * pow(R, n - (c * m - b - 1) / a); }

该文件通过#define MATRIX在「矩阵方式」（3×3 矩阵的Matrix<N>类型）与「贡献合并方式」（Info{x,y,s}类型，乘法为 $(x_1+x_2,y_1+y_2,s_1+s_2+x_2y_1)$）之间切换，两者均已在 Library Checker 验证通过。利用此模板，模板题 Library Checker - Sum of Floor of Linear 的实现只需分别给出 $U=(0,1,0)$、$R=(1,0,0)$（贡献合并版）或矩阵版 $U$、$R$，再调用euclid并取出结果的第三个分量。

例题四：【模板】类欧几里得算法（Luogu P5170，万能欧几里得版）

为应用万能欧几里得算法模板，首先将 $i=0$ 的项提出来单独考虑。剩余部分可看作对参数为 $(a,b,c,n)$ 的线段分别计算 $\sum y,\sum xy,\sum y^2$。与正文一致，有两种将操作序列转换为幺半群元素的方式。

矩阵运算。状态向量定义为 $(1,x,y,xy,y^2,\sum y,\sum xy,\sum y^2)$，初始状态为 $(1,0,0,0,0,0,0,0)$，两个操作分别为

$$ U = \begin{pmatrix} 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \ 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \ 0 & 0 & 1 & 0 & 2 & 0 & 0 & 0 \ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix},~ R = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \ 0 & 0 & 1 & 1 & 0 & 1 & 1 & 0 \ 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}. $$

最终答案为初始状态右乘这些操作矩阵的乘积得到的向量末尾三个分量。这一做法常数巨大，并不能通过本题，给出细节仅是为了辅助理解。

贡献合并。一段操作序列的贡献定义为 $(x,y,\sum y,\sum xy,\sum y^2)$，两个操作分别为

$$ U = (0,1,0,0,0),~ R = (1,0,0,0,0). $$

贡献合并时：

$$ \begin{aligned} \sum_{S_1+S_2} y &= \sum_{S_1}y + \sum_{S_2}(y+y_1) = \sum_{S_1}y + \sum_{S_2}y + x_2y_1,\ \sum_{S_1+S_2} xy &= \sum_{S_1}xy + \sum_{S_2}(x+x_1)(y+y_1) \ &= \sum_{S_1}xy + \sum_{S_2}xy + x_1\sum_{S_2}y + y_1\sum_{S_2}x + x_1y_1\sum_{S_2}1\ &= \sum_{S_1}xy + \sum_{S_2}xy + x_1\sum_{S_2}y + \dfrac{1}{2}x_2(x_2+1)y_1 + x_1x_2y_1,\ \sum_{S_1+S_2}y^2 &= \sum_{S_1}y^2 + \sum_{S_2}(y+y_1)^2 \ &= \sum_{S_1}y^2 + \sum_{S_2}y^2 + 2y_1\sum_{S_2}y + y_1^2\sum_{S_2}1 \ &= \sum_{S_1}y^2 + \sum_{S_2}y^2 + 2y_1\sum_{S_2}y + x_2y_1^2. \end{aligned} $$

这说明应将操作的乘法定义为

$$ \begin{aligned} &(x_1,y_1,s_1,t_1,u_1)\cdot(x_2,y_2,s_2,t_2,u_2)\ &= (x_1+x_2,y_1+y_2,s_1+s_2+x_2y_1,\ &\qquad t_1+t_2+x_1s_2+(1/2)x_2(x_2+1)y_1+x_1x_2y_1,\ &\qquad u_1+u_2+2y_1s_2+x_2y_1^2). \end{aligned} $$

虽然直接验证较为繁琐，但上述贡献向量在该乘法下确实构成幺半群，单位元为 $(0,0,0,0,0)$。

一般情形。有

$$ \begin{aligned} \sum_{S_1+S_2}x^ry^s &= \sum_{S_1}x^ry^s + \sum_{S_2}(x+x_1)^r(y+y_1)^s \ &= \sum_{S_1}x^ry^s + \sum_{i=0}^r\sum_{j=0}^s\binom{r}{i}\binom{s}{j}x_1^{r-i}y_1^{s-j}\sum_{S_2}x^iy^j. \end{aligned} $$

只要维护好所有更低幂次的贡献，就可以计算一般情形的和式。

仓库中的 euclidean-5.cpp 实现了本解法：Info结构维护五元组 $(x,y,s,t,u)$，在模 $998244353$ 下实现上述乘法（对应代码中的tmp = (rhs.x * (rhs.x + 1) / 2 + x * rhs.x) % M; res.t = ...; res.u = ...），并在主函数中用b / c单独补回 $i=0$ 项（f = res.s + b/c、h = res.u + (b/c)^2），输出顺序同样为f h g。

例题五：【清华集训 2014】Sum（Luogu P5172，万能欧几里得版）

单独处理 $r$ 为完全平方数的情形（与前文完全一致，从略），仅考虑 $r$ 非完全平方数的情形。

本题应用万能欧几里得算法的方式有很多。例如，可以为每个操作定义一个线性变换：

$$ U(x) = -x,~ R(x) = x + 1, $$

操作的乘法定义为线性变换的复合，最终答案就是操作序列对应的变换的复合函数在 $x=0$ 处的值。

还可以为每段操作序列定义贡献为 $((-1)^y,\sum(-1)^y)$，两个操作分别取

$$ U = (0,-1),~ R = (1,1), $$

贡献合并定义为

$$ (u_1,v_1)\cdot(u_2,v_2) = (u_1u_2,v_1+u_1v_2), $$

容易验证该乘法下所有操作构成幺半群，单位元为 $(0,1)$，最终答案是所有元素乘积的第二个分量。

这两种方法是一致的：如果将线性变换写作 $f(x)=u+vx$，那么线性变换复合对应的系数变化恰恰就是上述操作的乘法——这两个幺半群是同构的。

本题中线段的参数为 $(k,n)$，其中 $k\in\mathbf R$ 为直线斜率。设操作序列对应的乘积为 $F(k,n,U,R)$，递归算法如下：

• 若 $k\ge 1$，每个 $R$ 前方都有至少 $\lfloor k\rfloor$ 个 $U$，所以

$$ F(k,n,U,R) = F(k-\lfloor k\rfloor,n,U,U^{\lfloor k\rfloor} R). $$

• 若 $k<1$，交换操作序列中的 $U$ 和 $R$，并舍去末尾的 $U$（即交换前的 $R$），所以

$$ F(k,n,U,R) = F(k^{-1},m,R,U)R^{n-\lfloor k^{-1}m\rfloor}. $$

算法中 $k$ 的迭代过程其实就是在求 $\sqrt{r}$ 的连分数展开，为此可以应用 PQa 算法，求连分数的过程和万能欧几里得算法迭代的过程可以同时进行。和类欧几里得算法一致，算法复杂度仍是 $O(\log n)$ 的。

仓库中的 euclidean-6.cpp 实现了这一解法：LinearTransform{u,v}类型以eval(x)=u+v*x求值，乘法对应复合；主循环中P,Q维护二次无理数连分数展开的状态量（PQa 算法核心），a=(P+sqr)/Q为每轮连分数项，pow(U,a)*R完成 $U^{\lfloor k\rfloor}R$ 的合并，n=m完成横坐标收缩，并在n归零后退出循环。

习题推荐

模板题：

• Library Checker - Sum of Floor of Linear
• Luogu P5170【模板】类欧几里得算法
• Luogu P5171 Earthquake
• Luogu P5172 [清华集训 2014] Sum
• Luogu P4132 [BJOI2012] 算不出的等式
• LOJ 138. 类欧几里得算法
• LOJ 6440. 万能欧几里得
• Luogu P5179 Fraction
• Codeforces 1182 F. Maximum Sine

应用题：

• Luogu P4433 [COCI 2009/2010 #1] ALADIN
• AtCoder Beginner Contest 372 G - Ax + By < C
• AtCoder Beginner Contest 313 G - Redistribution of Piles
• AtCoder Beginner Contest 283 Ex - Popcount Sum
• Codeforces 1098 E. Fedya the Potter
• Codeforces 868 G. El Toll Caves

小结：两类算法的定位

类欧几里得算法与万能欧几里得算法共享同一个数学内核——分数的递归结构。前者以「取模 + 交换坐标」两步直接约化求和式，适合 $f,g,h$ 这类可写成（带权）直线下整点计数的和式，推导直接但每换一种和式都要重新推导；后者将操作序列抽象为幺半群元素，通过「提取整段 $U$ / 合并 $U^{\lfloor a/c\rfloor}R$ / 翻转坐标并交换 $U,R$」三种约化手段给出统一递归框架，只需替换幺半群类型T即可覆盖 $\sum y,\sum xy,\sum y^2$ 乃至一般 $\sum x^ry^s$ 的任意组合，甚至能处理斜率为无理数（二次无理数）的情形。在 OI 与 ICPC 实战中，掌握本文 euclidean.md 中的代数推导与几何直观，并能在 euclidean-4.cpp 的模板基础上定制T，即可覆盖上述模板题与应用题的绝大多数场景。

【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. （某大型游戏线上攻略，内含炫酷算术魔法）项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki