连分数
连分数
连分数 是实数作为有理数的特定收敛序列的表示。它们在竞争性编程(competitive programming)中很有用,因为它们易于计算,并且可以有效地用于在分母不超过给定值的所有数字中,找到基础实数(underlying real number)的最佳可能有理近似(best possible rational approximation)。
除此之外,连分数与欧几里得算法密切相关,这使得它们在一系列数论问题中非常有用。
定义
连分数是一种记号。例如,长为
只是为了形式上简洁,才记成等号左边的样子。这里的四个变元可以任意取值。
连分数各变元的下标从
简单连分数
可以证明,任何有理数都可以精确地以两种方式表示为连分数:
此外,对于
一旦深入研究了连分数构造的细节,这背后的原因就会很清楚。
定义
设
称为有理数
设
对于有限连分数,全体结尾为
简单连分数:连分数从第
简单连分数的值,一定大于偶数的渐进分数,一定小于奇数的渐进分数。无限简单连分数一定收敛。
仿照一般分数的概念,第
无限连分数
如果分式无限地写下去,有无限个变元,就得到无限连分数。无限连分数收敛等价于渐进分数收敛。
有定理:
无限连分数,如果各变元均大于等于
因为只要各变元为正,无限连分数的偶渐进分数单调递增(都比它小),奇渐进分数单调递减(都比它大)。而在均大于等于
显然可以看到,连分数关于下标为偶数的变元单调递增,关于下标为奇数的变元单调递减。这无论它有限或无限都成立。
定义
设
称为无理数
注意,对于
另一个重要的观察结果是,当
渐进分数
定义
在上面的定义中,有理数
相应地,单个
考虑
其中
请注意,在这种特定情况下,找到
定义
设
通常将大于
余项和部分商
定义
作为渐进分数的补充,定义 余项(完全商,complete quotients)为
相应地,将单个
相应地,取整得到的
对于各项均为正数的连分数,所有的余项也都是正数。
根据以上定义,可以得出
将
特别地,
这意味着可以从
序列
因此,对于任何无理数
求解简单连分数表示
在代码片段中主要假设有限的连分数。
如果要求
- 取倒数,得到的余项大于
。 - 取整得到整数部分为部分商,小数部分在
到 之间。 - 对小数部分重复上述操作。
这样就得到了相应的表示。
从
从这个表达式中,下一个完全商
对于
因此,
由此,
=== "C++"
```cpp
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;
}
```
=== "Python"
```py
def fraction(p, q):
a = []
while q:
a.append(p // q)
p, q = q, p % q
return a
```
如果规定第
如果两个无限简单连分数的值相等,必然逐项相等。
如果两个有限简单连分数的值相等,不仅要逐项相等,而且必然项数也相同。
无限简单连分数不能与有限简单连分数值相等。有理数与有限简单连分数具有一一对应关系,因此无限简单连分数全都是无理数。
性质
为了给连分数的进一步研究提供一些动力,现在给出一些性质。
递推
对于渐进分数
其中
渐进分数分子和分母具有完全相同的递推关系:
这里和 Farey 数列的递推关系很像。
形式上记初项:
只是形式上成立。第
可以注意到,
反序定理
如果
如果
对递推关系稍加改造,有:
又利用初值,即可证明反序定理。
渐进分数的差分
计算相邻两项渐进分数的差,需要通分。通分后的分子代入递推关系:
代入初值就有渐进分数的差分:
可以观察到,式
渐进分数的递推关系很像行列式的列变换。行列式一列加到另一列上不改变它的值,两列交换则反号。
根据递推式,如果连分数各项均为整数,则渐进分数分子分母总是互素。
对于有理数的简单连分数展开,常用渐进分数差分的等式,求解一次线性不定方程(参见 扩展欧几里得算法):
因为 a 与 b 互素,
倒数定理
由于实数与简单连分数一一对应,称实数的简单连分数的渐进分数,就是实数的渐进分数。于是就有倒数定理:
对于大于 1 的实数 x,x 的渐进分数的倒数恰好是
于是根据新的初值与递推就能发现倒数关系成立。
最佳逼近
让
那么
因此允许通过检查
下面会对这些性质建立一些直觉并做出进一步解释。
连行列式
继续看前面定义的渐近分数。对于
渐近分数是连分数的核心概念,因此研究它们的性质很重要。
对于数字
其中
因此,
为了一致性,定义了两个额外的渐近分数
渐进分数
根据渐进分数的定义,
由此得出
最初,
为了保持一致性,可以方便地定义
从数值分析可知,任意三对角矩阵的行列式
可以递归地计算为
这个多项式也被称为连行列式(continuant),由于其与连分数的密切关系。如果主对角线上的顺序颠倒,则连行列式不会改变。这产生了一个计算公式:
实现
把渐进分数计算为一对序列
=== "C++"
```cpp
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);
}
```
=== "Python"
```py
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 总小于其奇数阶渐近分数,大于其偶数阶渐近分数。
对于数字
特别地,这意味着
并且
由此可以得出结论
后一种不等式是由于
为了估计
将分子中的
因此,
因此
这产生了
根据递归关系,
总是定义明确的,因为基础系列总是收敛的。值得注意的是,剩余系列
由于
从这张图中可以看到
因此,
例题 扩展欧几里得
您将获得
虽然这个问题通常是用扩展欧几里得算法解决的,但有一个简单而直接的连分数的解决方案。
设
其中
```py
# 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]
```
几何解释
考虑线
奇数渐进分数
外壳上的所有整数顶点都作为
对于整数
在下图中可以看到
对于渐进分数
设
利用外积
最后一个等式是由于
注意到
正如已经注意到的,
在向量形式中,它重写为
这意味着
得出最终公式
例题 鼻子拉伸算法
每次将
因此,
换句话说,
在上面的图片中,
当不可能在不跨越
此过程生成接近直线的指数较长的向量。
对于这一特性,Boris Delaunay 将生成结果收敛向量的过程称为 鼻子拉伸算法(Nose stretching algorithm)。
如果观察在点
结合 Pick 定理,这意味着三角形内部没有严格的格点,其边界上的唯一格点是
这反过来意味着,具有奇数系数的
定义
这些多边形也被称为 克莱因多边形(Klein polygons),以费利克斯·克莱因(Felix Klein)的名字命名,他首次提出了对连续分数的几何解释。
例题
既然已经介绍了最重要的事实和概念,那么是时候深入研究具体的例题了。
线下凸包
找到格点
如果我们考虑无界集合
然而,在附加约束
设
然而,
为了到达外壳中的下一个格点,应该到达点
设
这样的点位于
也就是说,对于某些奇数
要找到这样的
当
并且
现在,可以将
```cpp
// 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);
}
```
=== "Python"
```py
# 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
您将得到整数
在这个问题中有
为了方便起见,通过替换
为了更一般地对待它,编写一个函数,该函数在
这个问题的核心解决方案思想基本上重复了前面的问题,但不是使用下中间分数来偏离直线,而是使用上中间分数来接近直线,而不跨越直线,也不违反
```py
# (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
计算
此和等于格点
在构造了
```cpp
// 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;
}
```
=== "Python"
```py
# 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
给定
如果您注意到
然而,将
```cpp
void solve(int p, int q, int N) {
cout << p * N * (N + 1) / 2 - q * sum_floor(fraction(p, q), N) << "\n";
}
```
=== "Python"
```py
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
给定
这是迄今为止技术上最麻烦的问题。
可以使用相同的方法来构造线
已经知道如何解决
现在应该注意到,一旦到达了离直线最近的点,就可以假设直线实际上通过了最近的点。因为在实际直线和稍微向下移动以通过最近点的直线之间,
也就是说,要在
```py
# 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
有一个有理数
这个问题等价于:查找
根据中国剩余定理,要求结果模化几个素数与要求其模化其乘积是相同的。因此,在不丧失一般性的情况下,假设知道余数模足够大的数
对于给定的余数
题面有
因此,问题归结为,给定
这实际上与找到
对于
由于
就连分数而言,这意味着
```py
# 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。
贡献者:@Great-designer@Menci@Chunibyo
本页面最近更新:2/3/2023, 12:00:00 AM,更新历史
发现错误?想一起完善? 在 GitHub 上编辑此页!
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用