类欧几里德算法
引入
类欧几里德算法是洪华敦在 2016 年冬令营营员交流中提出的内容。它常用于解决形如
⌊𝑎𝑖+𝑏𝑐⌋
结构的数列(下标为 𝑖
)的求和问题。它的主要想法是,利用分数自身的递归结构,将问题转化为更小规模的问题,递归求解。因为分数的递归结构和 欧几里得算法 存在直接的 联系,因此,这一求和方法也称为类欧几里得算法。
因为 连分数 和 Stern–Brocot 树 等方法同样刻画了分数的递归结构,所以利用类欧几里得算法可以解决的问题,通常也可以用这些方法解决。与这些方法相比,类欧几里得算法通常更容易理解,它的实现也更为简明。
类欧几里得算法
最简单的例子,就是求和问题:
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=0⌊𝑎𝑖+𝑏𝑐⌋,
其中,𝑎,𝑏,𝑐,𝑛
都是正整数。
代数解法
首先,将 𝑎,𝑏
对 𝑐
取模,可以简化问题,将问题转化为 0 ≤𝑎,𝑏 <𝑐
的情形:
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=0⌊𝑎𝑖+𝑏𝑐⌋=𝑛∑𝑖=0⌊(⌊𝑎𝑐⌋𝑐+(𝑎mod𝑐))𝑖+(⌊𝑏𝑐⌋𝑐+(𝑏mod𝑐))𝑐⌋=𝑛∑𝑖=0(⌊𝑎𝑐⌋𝑖+⌊𝑏𝑐⌋+⌊(𝑎mod𝑐)𝑖+(𝑏mod𝑐)𝑐⌋)=𝑛(𝑛+1)2⌊𝑎𝑐⌋+(𝑛+1)⌊𝑏𝑐⌋+𝑓(𝑎mod𝑐,𝑏mod𝑐,𝑐,𝑛).
现在,考虑转化后的问题。令
𝑚=⌊𝑎𝑛+𝑏𝑐⌋.
那么,原问题可以写作二次求和式:
𝑛∑𝑖=0⌊𝑎𝑖+𝑏𝑐⌋=𝑛∑𝑖=0𝑚−1∑𝑗=0[𝑗<⌊𝑎𝑖+𝑏𝑐⌋].
交换求和次序,这需要对于每个 𝑗
计算满足条件的 𝑖
的范围。为此,将条件变形:
𝑗<⌊𝑎𝑖+𝑏𝑐⌋=⌈𝑎𝑖+𝑏+1𝑐⌉−1⟺𝑗+1<⌈𝑎𝑖+𝑏+1𝑐⌉⟺𝑗+1<𝑎𝑖+𝑏+1𝑐⟺𝑐𝑗+𝑐−𝑏−1𝑎<𝑖⟺⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋<𝑖.
变形过程中多次利用了取整函数的性质。代入变形后的条件,原式可以写作:
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑚−1∑𝑗=0𝑛∑𝑖=0[𝑖>⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋]=𝑚−1∑𝑗=0(𝑛−⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋)=𝑛𝑚−𝑓(𝑐,𝑐−𝑏−1,𝑎,𝑚−1).
令 (𝑎′,𝑏′,𝑐′,𝑛′) =(𝑐,𝑐 −𝑏 −1,𝑎,𝑚 −1)
,这就又回到了前面讨论过的 𝑎′ >𝑐′
的情形。
将这两步转化结合在一起,可以发现在过程中,(𝑎,𝑐)
不断地取模后交换位置,直到 𝑎 =0
。这就类似于对 (𝑎,𝑐)
进行辗转相除,这也是类欧几里德算法的得名。它的时间复杂度是 𝑂(logmin{𝑎,𝑐})
的。
在计算过程中,可能会出现 𝑚 =0
的情形,此时内层递归会出现 𝑛 = −1
。这并不影响最终的结果。但是,如果要求出现 𝑚 =0
时,直接终止算法,那么算法的时间复杂度可以改良为 𝑂(logmin{𝑎,𝑐,𝑛})
的。
对复杂度的解释
利用该算法和欧几里得算法的相似性,很容易说明它的时间复杂度是 𝑂(logmin{𝑎,𝑐})
的。因此,只需要说明,如果在 𝑚 =0
时终止算法,那么它的时间复杂度也是 𝑂(log𝑛)
的。
令 𝑚 =⌊(𝑎𝑛 +𝑏)/𝑐⌋
,并记 𝑆 =𝑚𝑛
,𝑘 =𝑚/𝑛
,它们分别相当于几何直观(见下一节)中点阵图的面积和直线的斜率。对于充分大的 𝑛
,近似有 𝑘 ≐𝑎/𝑐
。
考察 𝑆
和 𝑘
在算法过程中的变化。第一步取模时,𝑛
保持不变,𝑘
近似由 𝑎/𝑐
变为 (𝑎mod𝑐)/𝑐
,相当于斜率由 𝑘
变为 𝑘 −⌊𝑘⌋
,而 𝑆
也近似变为原来的 (𝑘 −⌊𝑘⌋)
倍。第二步交换横纵坐标时,𝑆
近似保持不变,𝑘
则变为它的倒数。因此,若设两步操作后,二元组 (𝑘,𝑆)
变为 (𝑘′,𝑆′)
,则有 𝑘′ =(𝑘 −⌊𝑘⌋)−1
且 𝑆′ =(𝑘 −⌊𝑘⌋)𝑆
。
因为 1 ≤⌊𝑘′⌋ ≤𝑘′ <⌊𝑘′⌋ +1
,所以,递归计算两轮后,乘积缩小的倍数最少为
(𝑘′−⌊𝑘′⌋)(𝑘−⌊𝑘⌋)=1−⌊𝑘′⌋𝑘′<1−⌊𝑘′⌋⌊𝑘′⌋+1=1⌊𝑘′⌋+1≤12.
因此,至多 𝑂(log𝑆)
轮,算法必然终止。因为从第二轮开始,每轮开始时的 𝑆
总是不超过上一轮取模结束后的 𝑆
,而后者大致为 𝑘𝑛2
且 𝑘 <1
,因而 𝑂(log𝑆) ⊆𝑂(log𝑛)
。这就得到了上述结论。
模板题的参考实现如下:
模板题实现(Library Checker - Sum of Floor of Linear)
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;
}
|
几何直观
这个算法还可以从几何的角度理解。类欧几里得算法可以解决的问题主要是直线下整点计数问题。
如下图最左部分所示,该求和式相当于求直线
𝑦=𝑎𝑥+𝑏𝑐
下方,𝑥
轴上方(不包括 𝑥
轴),且横坐标位于 [0,𝑛]
之间的格点数目。

首先,移除斜率和截距中的整数部分。这一步相当于将上图中间部分的蓝点数量单独计算出来。当斜率和截距都是整数时,蓝点一定构成一个梯形阵列,也就是说,不同纵列的格点形成等差数列,因而这些点的数量是容易计算的。将这些点移除后,剩余的格点和上图最右部分的红点数量一致。问题就转化成了斜率和截距都小于一的情形。因为梯形的高为 𝑛 +1
且两个底边长度分别为 ⌊𝑏/𝑐⌋
和 (⌊𝑎/𝑐⌋𝑛 +⌊𝑏/𝑐⌋)
,所以,利用梯形面积公式,这一步骤可以归纳为算式
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑓(𝑎mod𝑐,𝑏mod𝑐,𝑐,𝑛)+12(𝑛+1)(⌊𝑏𝑐⌋+(⌊𝑎𝑐⌋𝑛+⌊𝑏𝑐⌋)).
然后,翻转横纵坐标轴。如下图最左部分所示,图中的红点和蓝点构成了一个横向长度为 𝑛
、纵向长度为 𝑚 =⌊(𝑎𝑛 +𝑏)/𝑐⌋
的矩形点阵。要计算红点的数量,只需要计算蓝点的数量,再用矩形点阵的数量减去蓝点的数量即可。翻转后,上图左半部分中的蓝点点阵就变成了某条直线下的红色点阵。而且,翻转后,斜率大于一,就又回到了上文已经处理过的情形。

关键在于如何计算新的红色点阵上方的直线的方程。将上图最左部分的横纵坐标轴翻转,就得到上图中间部分。翻转后的红色点阵上方的直线(中间部分的实线),并非对应翻转前的直线(最左部分的实线),而是翻转前的直线向左上平移一点点的结果(最左部分的虚线)。这是因为,如果直接将直线(最左部分的实线)翻转,将得到中间部分的虚线,但是按照定义,它下方的格点包括恰好落在直线上的格点,这就会导致直线上的格点重复计数。为了避免这一点,需要将翻转直线 𝑦 =(𝑎𝑥 +𝑏)/𝑐
后得到的直线 𝑦 =(𝑐𝑥 −𝑏)/𝑎
向下平移一点点,得到直线 𝑦 =(𝑐𝑥 −𝑏 −1)/𝑎
,这样它下方的点阵才恰为翻转前的蓝色点阵。
还有另一处细节需要处理。上图中间部分的直线的截距是负数,这意味着还没有回到之前的初始情形。要让截距恢复为非负数,只需要将直线(中间部分的实线)向左平移一个单位。这样做不会漏掉任何格点,因为翻转前的蓝色点阵中没有纵坐标为零的点,翻转后也就不存在横坐标为零的点。最后,直线方程就变为 𝑦 =(𝑐𝑥 +𝑐 −𝑏 −1)/𝑎
;同时,点阵的横坐标的上界也从 𝑚
变成了 𝑚 −1
。这一步骤可以归纳为算式
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑚𝑛−𝑓(𝑐,𝑐−𝑏−1,𝑎,𝑚−1).
这种递归的算法行得通,主要有两个原因:
- 第一,直线的斜率不断地先取小数部分再取倒数,这等价于计算直线斜率 𝑘 =𝑎/𝑐
的 连分数展开。因为有理分数的连分数展开的长度是 𝑂(logmin{𝑎,𝑐})
的,所以这一过程一定在 𝑂(logmin{𝑎,𝑐})
步后终止;
- 第二,因为每次翻转坐标轴的时候,直线斜率都是小于一的,因此,直觉上应该有 𝑚 <𝑛
,也就是说,经过这样一轮迭代后,横坐标的范围一直是在缩小的。前文的复杂度计算中通过严格的分析说明,每两轮迭代后,𝑛
至多为原来的一半,故而这一过程一定在 𝑂(log𝑛)
步后终止。
这也是斜率为有理数时的类欧几里得算法的复杂度是 𝑂(logmin{𝑎,𝑐,𝑛})
的原因。
利用类似的几何直观,可以将类欧几里得算法推广到斜率为无理数的情形,具体分析请参考后文的例题。
例题
【模板】类欧几里得算法
多组询问。给定正整数 𝑎,𝑏,𝑐,𝑛
,求
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=0⌊𝑎𝑖+𝑏𝑐⌋,𝑔(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=0𝑖⌊𝑎𝑖+𝑏𝑐⌋,ℎ(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=0⌊𝑎𝑖+𝑏𝑐⌋2.
解答一
类似于 𝑓
的推导,可以得到 𝑔,ℎ
的递归表达式。
首先,利用取模,将问题转化为 0 ≤𝑎,𝑏 <𝑐
的情形:
𝑔(𝑎,𝑏,𝑐,𝑛)=𝑔(𝑎mod𝑐,𝑏mod𝑐,𝑐,𝑛)+⌊𝑎𝑐⌋𝑛(𝑛+1)(2𝑛+1)6+⌊𝑏𝑐⌋𝑛(𝑛+1)2,ℎ(𝑎,𝑏,𝑐,𝑛)=ℎ(𝑎mod𝑐,𝑏mod𝑐,𝑐,𝑛)+2⌊𝑏𝑐⌋𝑓(𝑎mod𝑐,𝑏mod𝑐,𝑐,𝑛)+2⌊𝑎𝑐⌋𝑔(𝑎mod𝑐,𝑏mod𝑐,𝑐,𝑛)+⌊𝑎𝑐⌋2𝑛(𝑛+1)(2𝑛+1)6+⌊𝑏𝑐⌋2(𝑛+1)+⌊𝑎𝑐⌋⌊𝑏𝑐⌋𝑛(𝑛+1).
然后,利用交换求和次序,可以进一步转化。同样地,令
𝑚=⌊𝑎𝑛+𝑏𝑐⌋.
那么,对于和式 𝑔
,有
𝑔(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=0𝑖⌊𝑎𝑖+𝑏𝑐⌋=𝑛∑𝑖=0𝑚−1∑𝑗=0𝑖[𝑗<⌊𝑎𝑖+𝑏𝑐⌋]=𝑚−1∑𝑗=0𝑛∑𝑖=0𝑖[𝑖>⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋]=𝑚−1∑𝑗=012(⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋+𝑛+1)(𝑛−⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋)=12𝑚𝑛(𝑛+1)−12𝑚−1∑𝑗=0⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋−12𝑚−1∑𝑗=0⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋2=12𝑚𝑛(𝑛+1)−12𝑓(𝑐,𝑐−𝑏−1,𝑎,𝑚−1)−12ℎ(𝑐,𝑐−𝑏−1,𝑎,𝑚−1).
对于和式 ℎ
,有
ℎ(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=0⌊𝑎𝑖+𝑏𝑐⌋2=𝑛∑𝑖=0𝑚−1∑𝑗=0(2𝑗+1)[𝑗<⌊𝑎𝑖+𝑏𝑐⌋]=𝑚−1∑𝑗=0𝑛∑𝑖=0(2𝑗+1)[𝑖>⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋]=𝑚−1∑𝑗=0(2𝑗+1)(𝑛−⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋)=𝑛𝑚2−𝑚−1∑𝑗=0⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋−2𝑚−1∑𝑗=0𝑗⌊𝑐𝑗+𝑐−𝑏−1𝑎⌋=𝑛𝑚2−𝑓(𝑐,𝑐−𝑏−1,𝑎,𝑚−1)−2𝑔(𝑐,𝑐−𝑏−1,𝑎,𝑚−1).
从几何直观的角度看,这些非线性的求和式相当于给区域中的每个点 (𝑖,𝑗)
都赋予了相应的权重 𝑤(𝑖,𝑗)
。除了这些权重之外,其余部分的计算过程是完全一致的。对于权重的选择,一般地,有
𝑛∑𝑖=0𝑖𝑟⌊𝑎𝑖+𝑏𝑐⌋𝑠=𝑛∑𝑖=0𝑚−1∑𝑗=0𝑖𝑟((𝑗+1)𝑠−𝑗𝑠)[𝑗<⌊𝑎𝑖+𝑏𝑐⌋].
本题的另一个特点是,𝑔
和 ℎ
在递归计算时,会相互交错。因此,需要将 (𝑓,𝑔,ℎ)
作为三元组同时递归。
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(−1)⌊𝑑√𝑟⌋.
解答一
如果 𝑟
是完全平方数,那么当 √𝑟
为偶数时,和式为 𝑛
;否则,和式依据 𝑛
奇偶性不同,在 0
和 −1
之间交替变化。下面考虑 𝑟
不是完全平方数的情形。
为了应用类欧几里得算法,首先将求和式转化为熟悉的形式:
𝑛∑𝑑=1(−1)⌊𝑑√𝑟⌋=𝑛∑𝑑=1(1−2(⌊𝑑√𝑟⌋mod2))=𝑛−2𝑛∑𝑑=1(⌊𝑑√𝑟⌋−2⌊⌊𝑑√𝑟⌋2⌋)=𝑛−2𝑛∑𝑑=1⌊𝑑√𝑟⌋+4𝑛∑𝑑=1⌊𝑑√𝑟2⌋=𝑛−2𝑓(𝑛,1,0,1)+4𝑓(𝑛,1,0,2)
其中的函数 𝑓
具有形式
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=1⌊𝑎√𝑟+𝑏𝑐𝑖⌋.
与正文中的算法不同的是,此处的斜率不再是有理数。设斜率
𝑘=𝑎√𝑟+𝑏𝑐.
同样分为两种情形讨论。如果 𝑘 ≥1
,那么
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=1⌊𝑘𝑖⌋=𝑛∑𝑖=1⌊(𝑘−⌊𝑘⌋)𝑖⌋+⌊𝑘⌋𝑛∑𝑖=1𝑖=⌊𝑘⌋𝑛(𝑛+1)2+𝑓(𝑎,𝑏−𝑐⌊𝑘⌋,𝑐,𝑛).
问题转化为斜率小于一的情形。如果 𝑘 <1
,那么设 𝑚 =⌊𝑛𝑘⌋
,有
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=1⌊𝑘𝑖⌋=𝑛∑𝑖=1𝑚∑𝑗=1[𝑗≤⌊𝑘𝑖⌋]=𝑚∑𝑗=1𝑛∑𝑖=1[𝑖>⌊𝑘−1𝑗⌋]=𝑛𝑚−𝑚∑𝑗=1𝑛∑𝑖=1[𝑖≤⌊𝑘−1𝑗⌋].
此处的推导中,交换 𝑖
和 𝑗
的条件比正文中的情形更为简单,是因为直线 𝑦 =𝑘𝑥
上没有除了原点之外的格点。关键在于交换后的求和式写成 𝑓(𝑎,𝑏,𝑐,𝑛)
的形式,这相当于要求 𝑎′,𝑏′,𝑐′
满足
𝑘−1=𝑎′√𝑟+𝑏′𝑐′.
这并不困难,只需要将分母有理化,就能得到
𝑘−1=𝑐𝑎√𝑟+𝑏=𝑐𝑎√𝑟−𝑐𝑏𝑎2𝑟−𝑏2.
因此,有
𝑎′=𝑐𝑎, 𝑏′=−𝑐𝑏, 𝑐′=𝑎2𝑟−𝑏2.
这说明
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛𝑚−𝑓(𝑐𝑎,−𝑐𝑏,𝑎2𝑟−𝑏2,𝑚).
为了避免整数溢出,需要每次都将 𝑎,𝑏,𝑐
同除以它们的最大公约数。因为这个计算过程和计算 𝑘
的连分数的过程完全一致,所以根据 连分数理论,只要保证 gcd(𝑎,𝑏,𝑐) =1
,它们在计算过程中必然在整型范围内。另外,尽管 (𝑎,𝑏,𝑐,𝑛)
不会溢出,但是在该题数据范围下,𝑓(𝑎,𝑏,𝑐,𝑛)
可能会超过 64
位整数的范围,自然溢出即可,无需额外处理,最后结果一定在 [ −𝑛,𝑛]
之间。
尽管斜率不会变为零,算法的复杂度仍然是 𝑂(log𝑛)
的,这一点从前文关于算法复杂度的论证容易看出。
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;
}
|
Fraction
给定正整数 𝑎,𝑏,𝑐,𝑑
,求所有满足 𝑎/𝑏 <𝑝/𝑞 <𝑐/𝑑
的最简分数 𝑝/𝑞
中 (𝑞,𝑝)
的字典序最小的那个。
解答
这道题目也是 Stern–Brocot 树 的经典应用,相关题解可以在 此处 找到。因为它只依赖于分数的递归结构,所以它同样可以利用类似欧几里得算法的方法求解,故而也可以视作类欧几里得算法的一个应用。
如果 𝑎/𝑏
和 𝑐/𝑑
之间(不含端点)存在至少一个自然数,可以直接取 (𝑞,𝑝) =(1,⌊𝑎/𝑏⌋ +1)
。否则,必然有
⌊𝑎𝑏⌋≤𝑎𝑏<𝑝𝑞<𝑐𝑑≤⌊𝑎𝑏⌋+1.
从这个不等式中可以看出,𝑝/𝑞
的整数部分可以确定为 ⌊𝑎/𝑏⌋
,直接消去该整数部分,然后整体取倒数,用于确定它的小数部分。这正是确定 𝑝/𝑞
的连分数的 基本方法。若最终的答案是 𝑝/𝑞
,那么算法的时间复杂度为 𝑂(logmin{𝑝,𝑞})
。
此处,有一个细节问题,即取倒数之后得到的字典序最小的分数,是否是取倒数之前的字典序最小的分数。换句话说,满足 𝑎/𝑏 <𝑝/𝑞 <𝑐/𝑑
的分数 𝑝/𝑞
中,字典序 (𝑞,𝑝)
最小的,是否也是字典序 (𝑝,𝑞)
最小的。假设不然,设 𝑝/𝑞
是字典序 (𝑞,𝑝)
最小的,但是 𝑟/𝑠 ≠𝑝/𝑞
是字典序 (𝑟,𝑠)
最小的。这必然有 𝑟 <𝑝
且 𝑞 <𝑠
。但是,这说明
𝑎𝑏<𝑟𝑠<𝑟𝑞<𝑝𝑞<𝑐𝑑.
因此,𝑟/𝑞
无论按照哪个字典序怎样都是严格更小于当前解的。这与所设条件矛盾。因此,上述算法是正确的。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 | #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;
}
|
万能欧几里得算法
上一节讨论的类欧几里得算法推导通常较为繁琐,而且能够解决的和式主要是可以转化为直线下(带权)整点计数问题的和式。本节讨论一种更为一般的方法,它进一步抽象了上述过程,从而可以解决更多的问题。因此,这一方法也称为万能欧几里得算法。它同样利用了分数的递归结构求解问题,但是与类欧几里得算法约化问题的思路稍有不同。
仍然考虑最经典的求和问题:
𝑓(𝑎,𝑏,𝑐,𝑛)=𝑛∑𝑖=1⌊𝑎𝑖+𝑏𝑐⌋,
其中,𝑎,𝑏,𝑐,𝑛
都是正整数。
问题转化
设参数为 (𝑎,𝑏,𝑐,𝑛)
的线段为
𝑦=𝑎𝑥+𝑏𝑐, 0<𝑥≤𝑛.
对于这条线段,可以按照如下方法定义一个由 𝑈
和 𝑅
组成的字符串 𝑆
,也称为 操作序列:
- 字符串恰有 𝑛
个 𝑅
和 𝑚 =⌊(𝑎𝑛 +𝑏)/𝑐⌋
个 𝑈
组成;
- 第 𝑖
个 𝑅
前方的 𝑈
的数量恰等于 ⌊(𝑎𝑖 +𝑏)/𝑐⌋
,其中,𝑖 =1,⋯,𝑛
。
从几何直观上看,这大致相当于从原点开始,每向右穿过一次竖向的网格线,就写下一个 𝑅
,每向上穿过一次横向的网格线,就写下一个 𝑈
。如下图所示:

当然,这样的定义还需要考量一系列特殊情形:
- 经过整点(即同时上穿和右穿)时,需要先写 𝑈
再写 𝑅
;
- 字符串开始时,除了在 (0,1]
区间内上穿网格线的次数外,还需要格外补充 ⌊𝑏/𝑐⌋
个 𝑈
;
- 字符串结束时,不能有格外的 𝑈
。
如果对于几何直观的描述有任何不明晰的地方,可以参考上述代数方法的定义辅助理解。几何直观的描述,有助于理解下文的算法过程。
万能欧几里得算法的基本思路,就是将操作序列中的 𝑈
和 𝑅
都视作某个 幺半群 内的元素,将整个操作序列视为幺半群内元素的乘积,而问题最终的答案与这个乘积有关。
比如,本题中,可以定义状态向量 𝑣 =(1,𝑦,∑𝑦)
,表示自原点开始,经历了若干次上穿和右穿网格线后,当前的状态。其中,第一个分量是常数,第二个分量是纵坐标 𝑦
,第三个分量是要求的和式。起始时,有 𝑣 =(1,0,0)
。每向上穿过一次网格线,纵坐标就累加一,即相当于将状态向量右乘以矩阵
𝑈=⎛⎜
⎜
⎜⎝110010001⎞⎟
⎟
⎟⎠.
每向右穿过一次网格线,和式就累加一次纵坐标,即相当于将状态向量右乘以矩阵
𝑅=⎛⎜
⎜
⎜⎝100011001⎞⎟
⎟
⎟⎠.
因此,最终的状态就是乘积 (1,0,0)𝑆
,其中,𝑆
理解为上述矩阵的乘积。所求的答案,就是最终状态的第三个分量。
除了将幺半群中的元素定义为矩阵以外,还可以将它们定义为一段操作序列对于最终结果的贡献,然后将操作的乘积定义为两段操作序列的贡献的合并。
本题中,可以定义每段操作序列的贡献为 (𝑥,𝑦,∑𝑦)
。为了严谨地解释这些记号,可以将这些分量都看作是操作序列的函数,也就是说,对于操作序列 𝑆
,它的贡献可以写作 (𝑥(𝑆),𝑦(𝑆),(∑𝑦)(𝑆))
。其中,𝑥(𝑆)
和 𝑦(𝑆)
分别对应着操作序列 𝑆
中 𝑅
和 𝑈
的数量,也就是该线段右穿和上穿网格线的次数。最后一项中的求和符号,一般地,有如下定义:对于操作序列上的函数 𝑓(𝑆)
,可以定义 (∑𝑓)(𝑆)
,或记作 ∑𝑆𝑓
,为下面的表达式:
∑𝑆𝑓:=∑{𝑓(𝑆[1,𝑟]):𝑆𝑟=𝑅}.
其中,𝑆𝑟
是 𝑆
中的第 𝑟
个字符,𝑆[1,𝑟]
是 𝑆
中前 𝑟
个字符组成的前缀。也就是说,这个求和记号,可以看作是对于操作序列 𝑆
中所有以 𝑅
结尾的前缀进行的求和。比如说,有
∑𝑆1=𝑥, ∑𝑆𝑥=12𝑥(𝑥+1).
再比如说,∑𝑦
就是操作序列中,每次右穿网格线时,之前上穿网格线的次数的累加。对于整段操作序列来说,𝑦
在所有以 𝑅
结尾的前缀处的值,就是在 𝑖 =1,⋯,𝑛
处的所有 ⌊(𝑎𝑖 +𝑏)/𝑐⌋
的值。因此,对于整段操作序列计算的 ∑𝑦
,就是本题最终要求的量。
初始时,有 𝑈 =(0,1,0)
,𝑅 =(1,0,0)
。进一步,可以将两个元素 (𝑥1,𝑦1,𝑠1)
和 (𝑥2,𝑦2,𝑠2)
的乘积定义为
(𝑥1,𝑦1,𝑠1)⋅(𝑥2,𝑦2,𝑠2)=(𝑥1+𝑥2,𝑦1+𝑦2,𝑠1+𝑠2+𝑥2𝑦1).
其中,最后一项贡献合并的结果可以通过如下计算得到:
∑𝑆1+𝑆2𝑦=∑𝑆1𝑦+∑𝑆2(𝑦+𝑦1)=∑𝑆1𝑦+∑𝑆2𝑦+𝑦1∑𝑆21=𝑠1+𝑠2+𝑥2𝑦1.
容易验证,这个乘法运算满足结合律,且幺元为 (0,0,0)
,所以这些元素在该乘法运算下构成幺半群。所求的答案,就是乘积的第三个分量。
这两种方法都可以得到正确的结果。但是,因为保留了较多的冗余信息,矩阵运算的常数较大,所以第二种方法在处理实际问题时更为实用。
算法过程
与类欧几里得算法整体约化不同,万能欧几里得算法约化问题的手段是将这些操作分批次地合并。记字符串对应的操作的乘积为
𝐹(𝑎,𝑏,𝑐,𝑛,𝑈,𝑅).
约化过程具体如下:
-
当 𝑏 ≥𝑐
时,操作序列的开始有 ⌊𝑏/𝑐⌋
个 𝑈
,直接计算它们的乘积,并将这些 𝑈
从操作序列中移除。此时,第 𝑖
个 𝑅
前方的 𝑈
的数量等于
⌊𝑎𝑖+𝑏𝑐⌋−⌊𝑏𝑐⌋=⌊𝑎𝑖+(𝑏mod𝑐)𝑐⌋.
因此,这相当于将线段参数由 (𝑎,𝑏,𝑐,𝑛)
变为 (𝑎,𝑏mod𝑐,𝑐,𝑛)
。所以,对于这种情形,有
𝐹(𝑎,𝑏,𝑐,𝑛,𝑈,𝑅)=𝑈⌊𝑏/𝑐⌋𝐹(𝑎,𝑏mod𝑐,𝑐,𝑛,𝑈,𝑅).
-
当 𝑎 ≥𝑐
时,操作序列中每个 𝑅
的前方都至少有 ⌊𝑎/𝑐⌋
个 𝑈
,可以将它们合并到 𝑅
上。也就是说,可以用 𝑈⌊𝑎/𝑐⌋𝑅
替代 𝑅
。合并后的字符串中,第 𝑖
个 𝑅
前方的 𝑈
的数量等于
⌊𝑎𝑖+𝑏𝑐⌋−⌊𝑎𝑐⌋𝑖=⌊(𝑎mod𝑐)𝑖+𝑏𝑐⌋.
因此,这相当于将线段参数由 (𝑎,𝑏,𝑐,𝑛)
变为 (𝑎mod𝑐,𝑏,𝑐,𝑛)
。所以,对于这种情形,有
𝐹(𝑎,𝑏,𝑐,𝑛,𝑈,𝑅)=𝐹(𝑎mod𝑐,𝑏,𝑐,𝑛,𝑈,𝑈⌊𝑎/𝑐⌋𝑅).
-
对于剩下的情形,需要翻转横纵坐标,这基本是在交换 𝑈
和 𝑅
,只是翻转后线段的参数需要仔细计算。结合操作序列的定义可知,需要确定系数 (𝑎′,𝑏′,𝑐′,𝑛′)
使得变换前的操作序列中,第 𝑗
个 𝑈
前方的 𝑅
的数量恰为 ⌊(𝑎′𝑗 +𝑏′)/𝑐′⌋
且总共有 𝑛′
个 𝑈
。根据定义可知,
𝑛′=⌊𝑎𝑛+𝑏𝑐⌋=𝑚,
而第 𝑗
个 𝑈
前方的 𝑅
的数量,就等于最大的 𝑖
使得
⌊𝑎𝑖+𝑏𝑐⌋<𝑗⟺𝑎𝑖+𝑏𝑐<𝑗⟺𝑖<𝑐𝑗−𝑏𝑎⟺𝑖<⌈𝑐𝑗−𝑏𝑎⌉=⌊𝑐𝑗−𝑏−1𝑎⌋+1.
因此,𝑖 =⌊(𝑐𝑗 −𝑏 −1)/𝑎⌋
。这一推导过程与前文类欧几里得算法的推导类似,同样利用了上下取整函数的性质。
有两处细节需要处理:
- 截距项 −(𝑏 +1)/𝑎
为负数。注意到,如果将线段向左平移一个单位,就可以让截距项恢复为非负数,因为总有 (𝑐 −𝑏 −1)/𝑎 ≥0
。因此,可以将交换前的第一段 𝑅⌊(𝑐−𝑏−1)/𝑎⌋𝑈
提取出来,只交换剩余操作序列中的 𝑈
和 𝑅
;
- 交换 𝑈
和 𝑅
后,结尾存在多余的 𝑈
。因此,交换 𝑈
和 𝑅
之前,需要首先将最后一段 𝑅
提取出来,只交换剩余操作序列中的 𝑈
和 𝑅
。这一段 𝑅
的数量为 𝑛 −⌊(𝑐𝑚 −𝑏 −1)/𝑎⌋
。
去掉头尾若干个字符后,第 𝑗
个 𝑈
前方的 𝑅
的数量变为:
⌊𝑐(𝑗+1)−𝑏−1𝑎⌋−⌊𝑐−𝑏−1𝑎⌋=⌊𝑐𝑗+(𝑐−𝑏−1)mod𝑎𝑎⌋.
回忆起,交换前的序列中 𝑈
的数量为 𝑚 =⌊(𝑎𝑛 +𝑏)/𝑐⌋
。而上述左移一个单位的操作,要求保证交换前至少存在一个 𝑈
,也就是 𝑚 >0
。利用这一条件,可以分为两种情形:
-
对于 𝑚 >0
的情形,处理了上面的两点后,交换完 𝑈
和 𝑅
的操作序列就是对应着参数为 (𝑐,(𝑐 −𝑏 −1)mod𝑎,𝑎,𝑚 −1)
的线段的合法序列。所以,有
𝐹(𝑎,𝑏,𝑐,𝑛,𝑈,𝑅)=𝑅⌊(𝑐−𝑏−1)/𝑎⌋𝑈𝐹(𝑐,(𝑐−𝑏−1)mod𝑎,𝑎,𝑚−1,𝑅,𝑈)𝑅𝑛−⌊(𝑐𝑚−𝑏−1)/𝑎⌋.
-
特别地,对于 𝑚 =0
的情形,交换前的操作序列中只包含 𝑛
个 𝑅
,无需交换,可以直接返回:
𝐹(𝑎,𝑏,𝑐,𝑛,𝑈,𝑅)=𝑅𝑛.
与类欧几里得算法不同,万能欧几里得算法的这一特殊情形需要单独处理,否则会因涉及负幂次而无法正确计算。
利用这些讨论,就可以将问题递归地解决。
假设幺半群内元素单次相乘的时间复杂度是 𝑂(1)
的。那么,如果计算过程中这些元素的幂次计算都使用 快速幂 进行,最终的算法复杂度就是 𝑂(logmax{𝑎,𝑐} +log(𝑏/𝑐))
的[^complexity]。
对复杂度的解释
对比(类)欧几里得算法,万能欧几里得算法只是多了求快速幂的步骤。其余的计算过程的复杂度和类欧几里得算法相仿,已经说明是 𝑂(logmin{𝑎,𝑐,𝑛})
的。现在,需要计算这些快速幂的总复杂度。
除了第一轮迭代,都有 $b
本页面最近更新:2025/6/1 02:56:59,更新历史
发现错误?想一起完善? 在 GitHub 上编辑此页!
本页面贡献者:sshwy, StudyingFather, Enter-tainer, H-J-Granger, Tiphereth-A, countercurrent-time, NachtgeistW, Early0v0, Ir1d, MegaOwIer, Xeonacid, AngelKitty, CCXXXI, cjsoft, diauweb, ezoixx130, FFjet, GekkaSaori, Henry-ZHR, Konano, LovelyBuggies, Makkiy, mgt, minghu6, P-Y-Y, PotassiumWings, SamZhangQingChuan, Suyun514, weiyong1024, alphagocc, c-forrest, cxm1024, GavinZhengOI, Gesrua, Great-designer, iamtwz, ksyx, kxccc, lychees, megakite, Peanut-Tang, r-value, SukkaW, TonyYin0418
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用