连分数

连分数

连分数 是实数作为有理数的特定收敛序列的表示。它们在竞争性编程（competitive programming）中很有用，因为它们易于计算，并且可以有效地用于在分母不超过给定值的所有数字中，找到基础实数（underlying real number）的最佳可能有理近似（best possible rational approximation）。

除此之外，连分数与欧几里得算法密切相关，这使得它们在一系列数论问题中非常有用。

定义

连分数是一种记号。例如，长为 \(4\) 的连分数：

只是为了形式上简洁，才记成等号左边的样子。这里的四个变元可以任意取值。

连分数各变元的下标从 \(0\) 开始。

简单连分数

可以证明，任何有理数都可以精确地以两种方式表示为连分数：

此外，对于 \(r=\frac{p}{q}\)，这种连分数的长度 \(k\) 估计为 \(k = O(\log \min(p, q))\)。

一旦深入研究了连分数构造的细节，这背后的原因就会很清楚。

定义

设 \(a_0, a_1, \dots, a_k \in \mathbb Z\) 和 \(a_1, a_2, \dots, a_k \geq 1\)。然后表达式

称为有理数 \(r\) 的 连分数表示，并简短地表示为 \(r=[a_0;a_1,a_2,\dots,a_k]\)。

对于有限连分数，全体结尾为 \(1\) 的有限连分数和全体结尾不为 \(1\) 的有限连分数一一对应，即同一个连分数有两种表示：

\([a_0,a_1,a_2,a_3]=[a_0,a_1,a_2,a_3-1,1]\)

简单连分数：连分数从第 \(1\) 项开始全都是正整数。如果有限，要求最后一项不为 \(1\)。（第 \(0\) 项可以任意）

简单连分数的值，一定大于偶数的渐进分数，一定小于奇数的渐进分数。无限简单连分数一定收敛。

仿照一般分数的概念，第 \(0\) 项是 \(0\) 的连分数称为「真分数」。显然如果这之后的所有变元都大于等于 \(1\)，那么得到的真分数一定落在 \(0\) 到 \(1\) 之间。

无限连分数

如果分式无限地写下去，有无限个变元，就得到无限连分数。无限连分数收敛等价于渐进分数收敛。

有定理：

无限连分数，如果各变元均大于等于 \(1\)，那么一定收敛。

因为只要各变元为正，无限连分数的偶渐进分数单调递增（都比它小），奇渐进分数单调递减（都比它大）。而在均大于等于 \(1\) 时，相邻（奇偶间）两个渐进分数之间距离可以给出估计式，趋于 \(0\)，因此收敛。

显然可以看到，连分数关于下标为偶数的变元单调递增，关于下标为奇数的变元单调递减。这无论它有限或无限都成立。

定义

设 \(a_0,a_1,a_2, \dots\) 为整数序列，使得 \(a_1, a_2, \dots \geq 1\)。设 \(r_k = [a_0; a_1, \dots, a_k]\)。然后表达式

称为无理数 \(r\) 的 连分数表示，并简短地表示为 \(r = [a_0;a_1,a_2,\dots]\)。

注意，对于 \(r=[a_0;a_1,\dots]\) 和整数 \(k\)，有 \(r+k = [a_0+k; a_1, \dots]\)。

另一个重要的观察结果是，当 \(a_0>0\) 时，\(\frac{1}{r}=[0;a_0, a_1, \dots]\)；当 \(a_1＝0\) 时，\(\frac{1}{r} = [a_1; a_2, \dots]\)。

渐进分数

定义

在上面的定义中，有理数 \(r_0, r_1, r_2, \dots\) 称为 \(r\) 的 渐进分数（convergents，意为「收敛」）。

相应地，单个 \(r_k = [a_0; a_1, \dots, a_k] = \frac{p_k}{q_k}\) 称为 \(r\) 的第 \(k\) 个渐进分数。

定义

设 \(r_k = [a_0; a_1, \dots, a_{k-1}, a_k]\)。对于 \(1 \leq t \leq a_k\)，\([a_0; a_1, \dots, a_{k-1}, t]\) 称为 中间分数（semiconvergents，「semi」意为「半」）。

通常将大于 \(r\) 的分数称为 上（upper）渐进分数或中间分数，将小于 \(r\) 者称为 下（lower）渐进分数或中间分数。

余项和部分商

定义

作为渐进分数的补充，定义 余项（完全商，complete quotients）为 \(s_k = [a_k; a_{k+1}, a_{k+2}, \dots]\)。

相应地，将单个 \(s_k\) 称为 \(r\) 的第 \(k\) 个完全商。

相应地，取整得到的 \(a_k\) 称为部分商。

对于各项均为正数的连分数，所有的余项也都是正数。

根据以上定义，可以得出 \(s_k \geq 1\) 代表 \(k \geq 1\) 的结论。

将 \([a_0; a_1, \dots, a_k]\) 视为一个形式代数表达式，并允许任意实数代替 \(a_i\)，得到

特别地，\(r = [s_0] = s_0\)。另一方面，可以将 \(s_k\) 表示为

这意味着可以从 \(s_k\) 计算 \(a_k = \lfloor s_k \rfloor\) 和 \(s_{k+1} = (s_k - a_k)^{-1}\)。

序列 \(a_0, a_1, \dots\) 定义良好，除非 \(s_k=a_k\)，这仅在 \(r\) 为有理数时发生。

因此，对于任何无理数 \(r\)，连分数表示都是唯一的。

求解简单连分数表示

在代码片段中主要假设有限的连分数。

如果要求 \((0,1)\) 区间内某个数的简单连分数表示（第 \(0\) 项为 \(0\)），只需：

• 取倒数，得到的余项大于 \(1\)。
• 取整得到整数部分为部分商，小数部分在 \(0\) 到 \(1\) 之间。
• 对小数部分重复上述操作。

这样就得到了相应的表示。

从 \(s_k\) 到 \(s_{k+1}\) 的转换如下

从这个表达式中，下一个完全商 \(s_{k+1}\) 如下

对于 \(s_k=\frac{p}{q}\)，这意味着

因此，\(r=\frac{p}{q}\) 的连分数表示的计算遵循 \(p\) 和 \(q\) 的欧几里得算法的步骤。

由此，\(\frac{p_k}{q_k} = [a_0; a_1, \dots, a_k]\) 的 \(\gcd(p_k, q_k) = 1\)。因此，渐进分数总是不可约的。

1
2
3
4
5
6
7
8

auto fraction(int p, int q) {
    vector<int> a;
    while(q) {
        a.push_back(p / q);
        tie(p, q) = make_pair(q, p % q);
    }
    return a;
}

1
2
3
4
5
6

def fraction(p, q):
    a = []
    while q:
        a.append(p // q)
        p, q = q, p % q
    return a

如果规定第 \(0\) 项是该数的取整，那么全体实数都有「唯一的简单连分数表示」。其中：

如果两个无限简单连分数的值相等，必然逐项相等。

如果两个有限简单连分数的值相等，不仅要逐项相等，而且必然项数也相同。

无限简单连分数不能与有限简单连分数值相等。有理数与有限简单连分数具有一一对应关系，因此无限简单连分数全都是无理数。

性质

为了给连分数的进一步研究提供一些动力，现在给出一些性质。

递推

渐进分数分子和分母具有完全相同的递推关系：

这里和 Farey 数列的递推关系很像。

形式上记初项：

只是形式上成立。第 \(-1\) 项渐进分数是 1/0，没有实际意义。

反序定理

如果 \(a_0\neq 0\)：

如果 \(a_0=0\)：

渐进分数的差分

计算相邻两项渐进分数的差，需要通分。通分后的分子代入递推关系：

代入初值就有渐进分数的差分：

根据递推式，如果连分数各项均为整数，则渐进分数分子分母总是互素。

对于有理数的简单连分数展开，常用渐进分数差分的等式，求解一次线性不定方程（参见 扩展欧几里得算法）：

因为 a 与 b 互素，\(\frac{a}{b}\) 就是最简的有理数，也就是它本身的最后一个渐进分数。那么，它的前一个渐进分数就是所求的解。

倒数定理

由于实数与简单连分数一一对应，称实数的简单连分数的渐进分数，就是实数的渐进分数。于是就有倒数定理：

对于大于 1 的实数 x，x 的渐进分数的倒数恰好是 \(\frac{1}{x}\) 的渐进分数。显然，该定理也应该对于 0 到 1 之间的实数 x 成立。

最佳逼近

因此允许通过检查 \(r\) 是否中间分数来找到其最佳有理逼近。

下面会对这些性质建立一些直觉并做出进一步解释。

连行列式

继续看前面定义的渐近分数。对于 \(r=[a_0, a_1, a_2, \dots]\)，其渐近分数为

渐近分数是连分数的核心概念，因此研究它们的性质很重要。

对于数字 \(r\)，其第 \(k\) 个渐近分数 \(r_k = \frac{p_k}{q_k}\) 可以计算为

其中 \(P_k(a_0,\dots,a_k)\) 是连行列式（continuant），定义为

因此，\(r_k\) 是 \(r_{k-1}\) 和 \(r_{k-2}\) 的加权中间值（mediant）。

为了一致性，定义了两个额外的渐近分数 \(r_{-1} = \frac{1}{0}\) 和 \(r_{-2} = \frac{0}{1}\)。

实现

把渐进分数计算为一对序列 \(p_{-2}, p_{-1}, p_0, p_1, \dots, p_k\) 和 \(q_{-2}, q_{-1}, q_0, q_1, \dots, q_k\)：

1
2
3
4
5
6
7
8
9

auto convergents(vector<int> a) {
    vector<int> p = {0, 1};
    vector<int> q = {1, 0};
    for(auto it: a) {
        p.push_back(p[p.size() - 1] * it + p[p.size() - 2]);
        q.push_back(q[q.size() - 1] * it + q[q.size() - 2]);
    }
    return make_pair(p, q);
}

1
2
3
4
5
6
7

def convergents(a):
    p = [0, 1]
    q = [1, 0]
    for it in a:
        p.append(p[-1]*it + p[-2])
        q.append(q[-1]*it + q[-2])
    return p, q

误差和余项的估计

实数 x 也可以写成：

最后一项渐近分数就是 x 本身。于是根据渐进分数的递推式，就有：

于是可以估计渐进分数的误差：

分别对 k 取奇数偶数就得到，x 总小于其奇数阶渐近分数，大于其偶数阶渐近分数。

对于数字 \(r\) 及其第 \(k\) 个渐进分数 \(r_k=\frac{p_k}{q_k}\)，以下公式成立：

特别地，这意味着

并且

由此可以得出结论

后一种不等式是由于 \(r_k\) 和 \(r_{k+1}\) 通常位于 \(r\) 的不同侧面，因此

例题 扩展欧几里得

您将获得 \(A, B, C \in \mathbb Z\)。查找 \(x, y \in \mathbb Z\)，使 \(Ax + By = C\).

1
2
3
4
5
6
7

# return (x, y) such that Ax+By=C
# assumes 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]

几何解释

对于渐进分数 \(r_k = \frac{p_k}{q_k}\)，设 \(\vec r_k = (q_k;p_k)\)。然后，以下重复出现：

设 \(\vec r = (1;r)\)。然后，每个向量 \((x;y)\) 对应于等于其斜率系数 \(\frac{y}{x}\) 的数字。

利用外积 \((x_1;y_1) \times (x_2;y_2) = x_1 y_2 - x_2 y_1\) 的概念，可以看出（参见下面的解释）

最后一个等式是由于 \(r_{k-1}\) 和 \(r_{k-2}\) 位于 \(r\) 的不同侧，因此 \(\vec r_{k-1}\) 和 \(\vec r_{k-2}\) 与 \(\vec r\) 的外积具有不同的符号。考虑到 \(a_k = \lfloor s_k \rfloor\)，\(\vec r_k\) 的公式如下

注意到 \(\vec r_k \times r = (q;p) \times (1;r) = qr - p\)，因此

例题 鼻子拉伸算法

每次将 \(\vec r_{k-1}\) 添加到向量 \(\vec p\) 时，\(\vec p \times \vec r\) 的值都会增加 \(\vec r_{k-1} \times \vec r\)。

因此，\(a_k=\lfloor s_k \rfloor\) 是 \(\vec r_{k-1}\) 向量的最大整数，可以将其添加到 \(\vec r_{k-2}\)，而无需更改与 \(\vec r\) 的外积的符号。

换句话说，\(a_k\) 是您可以将 \(\vec r_{k-1}\) 添加到 \(\vec r_{k-2}\) 的最大整数次数，而无需跨越 \(\vec r\) 定义的线：

在上面的图片中，\(\vec r_2 = (4;3)\) 是通过将 \(\vec r_1 = (1;1)\) 重复添加到 \(\vec r_0 = (1;0)\) 而获得的。

当不可能在不跨越 \(y=rx\) 线的情况下将 \(\vec r_1\) 进一步添加到 \(\vec r_0\) 时，转到另一侧，重复将 \(\vec r_2\) 添加到 \(\vec r_1\) 以获得 \(\vec r_3 = (9;7)\)。

此过程生成接近直线的指数较长的向量。

对于这一特性，Boris Delaunay 将生成结果收敛向量的过程称为 鼻子拉伸算法（Nose stretching algorithm）。

如果观察在点 \(\vec r_{k-2}\)、\(\vec r_{k}\) 和 \(\vec 0\) 上绘制的三角形，会注意到它的加倍面积是

结合 Pick 定理，这意味着三角形内部没有严格的格点，其边界上的唯一格点是 \(\vec 0\) 和 \(\vec r_{k-2} + t \cdot \vec r_{k-1}\)，对于所有整数 \(t\)，使得 \(0 \leq t \leq a_k\)。当连接所有可能的 \(k\) 时，这意味着在由偶数索引和奇数索引收敛向量形成的多边形之间的空间中没有整数点。

这反过来意味着，具有奇数系数的 \(\vec r_k\) 形成了线 \(y=rx\) 上方 \(x \geq 0\) 的格点凸包，而具有偶数系数的 \(\vec r_k\) 形成线 \(y=rx\) 下方 \(x > 0\) 的格点凸包。

定义

这些多边形也被称为 克莱因多边形（Klein polygons），以费利克斯·克莱因（Felix Klein）的名字命名，他首次提出了对连续分数的几何解释。

例题

既然已经介绍了最重要的事实和概念，那么是时候深入研究具体的例题了。

线下凸包

找到格点 \((x;y)\) 的凸包，使得 \(r=[a_0;a_1,\dots,a_k]=\frac{p_k}{q_k}\) 的 \(0 \leq x \leq N\) 和 \(0 \leq y \leq rx\)。

 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

// returns [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(auto a, int N) {
    auto [p, q] = convergents(a);
    int t = N / q.back();
    vector ah = {t};
    vector ph = {0, t*p.back()};
    vector qh = {0, t*q.back()};

    for(int i = q.size() - 1; i >= 0; 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

# returns [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

Timus - Crime and Punishment

您将得到整数 \(A\)、\(B\) 和 \(N\)。查找 \(x \geq 0\) 和 \(y \geq 0\)，使 \(Ax + By \leq N\) 和 \(Ax + By\) 达到最大值。

 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

# (x, y) such that y = (A*x+B) // C,
# Cy - Ax is max and 0 <= x <= N.
def closest(A, B, C, N):
    # y <= (A*x + B)/C <=> diff(x, y) <= B
    def diff(x, y):
        return C*y-A*x
    a = fraction(A, C)
    p, q = convergents(a)
    ph = [B // C]
    qh = [0]
    for i in range(2, len(q) - 1):
        if i % 2 == 0:
            while diff(qh[-1] + q[i+1], ph[-1] + p[i+1]) <= B:
                t = 1 + (diff(qh[-1] + q[i-1], ph[-1] + p[i-1]) - B - 1) // abs(diff(q[i], p[i]))
                dp = p[i-1] + t*p[i]
                dq = q[i-1] + t*q[i]
                k = (N - qh[-1]) // dq
                if k == 0:
                    return qh[-1], ph[-1]
                if diff(dq, dp) != 0:
                    k = min(k, (B - diff(qh[-1], ph[-1])) // diff(dq, dp))
                qh.append(qh[-1] + k*dq)
                ph.append(ph[-1] + k*dp)
    return qh[-1], ph[-1]

def solve(A, B, N):
    x, y = closest(A, N % A, B, N // A)
    return N // A - x, y

June Challenge 2017 - Euler Sum

计算 \(\sum\limits_{x=1}^N \lfloor \mathrm{e}x \rfloor\)，其中 \(\mathrm{e} = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, \dots, 1, 2n, 1, \dots]\) 是自然对数的底，\(N \leq 10^{4000}\)。

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

// sum floor(k * x) for k in [1, N] and x = [a0; a1, a2, ...]
int sum_floor(auto a, int N) {
    N++;
    auto [ah, ph, qh] = hull(a, N);

    // The number of lattice points within a vertical right trapezoid
    // on points (0; 0) - (0; y1) - (dx; y2) - (dx; 0) that has
    // a+1 integer points on the segment (0; y1) - (dx; y2).
    auto picks = [](int y1, int y2, int dx, int a) {
        int b = y1 + y2 + a + dx;
        int A = (y1 + y2) * dx;
        return (A - b + 2) / 2 + b - (y2 + 1);
    };

    int ans = 0;
    for(size_t i = 1; i < qh.size(); i++) {
        ans += picks(ph[i - 1], ph[i], qh[i] - qh[i - 1], ah[i - 1]);
    }
    return ans - N;
}

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

# sum floor(k * x) for k in [1, N] and x = [a0; a1, a2, ...]
def sum_floor(a, N):
    N += 1
    ah, ph, qh = hull(a, N)

    # The number of lattice points within a vertical right trapezoid
    # on points (0; 0) - (0; y1) - (dx; y2) - (dx; 0) that has
    # a+1 integer points on the segment (0; y1) - (dx; y2).
    def picks(y1, y2, dx, a):
        b = y1 + y2 + a + dx
        A = (y1 + y2) * dx
        return (A - b + 2) // 2 + b - (y2 + 1)

    ans = 0
    for i in range(1, len(qh)):
        ans += picks(ph[i-1], ph[i], qh[i]-qh[i-1], ah[i-1])
    return ans - N

NAIPC 2019 - It's a Mod, Mod, Mod, Mod World

给定 \(p\)、\(q\) 和 \(n\)，计算 \(\sum\limits_{i=1}^n [p \cdot i \bmod q]\)。

1
2
3

void solve(int p, int q, int N) {
    cout << p * N * (N + 1) / 2 - q * sum_floor(fraction(p, q), N) << "\n";
}

1
2

def solve(p, q, N):
    return p * N * (N + 1) // 2 - q * sum_floor(fraction(p, q), N)

Library Checker - Sum of Floor of Linear

给定 \(N\)、\(M\)、\(A\) 和 \(B\)，计算 \(\sum\limits_{i=0}^{N-1} \lfloor \frac{A \cdot i + B}{M} \rfloor\)。

 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

# hull of lattice (x, y) such that C*y <= A*x+B
def hull(A, B, C, N):
    def diff(x, y):
        return C*y-A*x
    a = fraction(A, C)
    p, q = convergents(a)
    ah = []
    ph = [B // C]
    qh = [0]

    def insert(dq, dp):
        k = (N - qh[-1]) // dq
        if diff(dq, dp) > 0:
            k = min(k, (B - diff(qh[-1], ph[-1])) // diff(dq, dp))
        ah.append(k)
        qh.append(qh[-1] + k*dq)
        ph.append(ph[-1] + k*dp)

    for i in range(1, len(q) - 1):
        if i % 2 == 0:
            while diff(qh[-1] + q[i+1], ph[-1] + p[i+1]) <= B:
                t = (B - diff(qh[-1] + q[i+1], ph[-1] + p[i+1])) // abs(diff(q[i], p[i]))
                dp = p[i+1] - t*p[i]
                dq = q[i+1] - t*q[i]
                if dq < 0 or qh[-1] + dq > N:
                    break
                insert(dq, dp)

    insert(q[-1], p[-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]
                insert(dq, dp)
    return ah, ph, qh

OKC 2 - From Modular to Rational

有一个有理数 \(\frac{p}{q}\)，即 \(1 \leq p, q \leq 10^9\)。您可以询问几个素数 \(m\) 的 \(p q^{-1}\) 模 \(m \sim 10^9\) 的值。恢复 \(\frac{p}{q}\)。

这个问题等价于：查找 \(1 \leq x \leq N\) 中，使 \(Ax \bmod M\) 最小的 \(x\)。

1
2
3
4
5
6
7
8

# find Q that minimizes Q*r mod m for 1 <= k <= n < m 
def mod_min(r, n, m):
    a = fraction(r, m)
    p, q = convergents(a)
    for i in range(2, len(q)):
        if i % 2 == 1 and (i + 1 == len(q) or q[i+1] > n):
            t = (n - q[i-1]) // q[i]
            return q[i-1] + t*q[i]

习题

• UVa OJ - Continued Fractions
• ProjectEuler+ #64: Odd period square roots
• Codeforces Round #184 (Div. 2) - Continued Fractions
• Codeforces Round #201 (Div. 1) - Doodle Jump
• Codeforces Round #325 (Div. 1) - Alice, Bob, Oranges and Apples
• POJ Founder Monthly Contest 2008.03.16 - A Modular Arithmetic Challenge
• 2019 Multi-University Training Contest 5 - fraction
• SnackDown 2019 Elimination Round - Election Bait

本页面主要译自博文 Continued fractions，版权协议为 CC-BY-SA 4.0。

本页面的全部内容在 小熊老师 - 莆田青少年编程俱乐部 0594codes.cn 协议之条款下提供，附加条款亦可能应用