让我们从一个定义开始,准确地确定在任何特定情况下我们正在寻找的分数:
定义。说一个分数a/b比另一个分数c/d更简单(两个分数都用最低的术语写成,
正分母)如果b <= d,abs(a) <= abs(c),至少
这两个不等式之一是严格的。
例如5/7 比6/7 简单,5/7 比5/8 简单,但2/5 和3/4 都不比另一个简单。 (我们这里没有总订购量。)
然后有了这个定义,就有一个非立即显而易见的定理,它保证我们正在寻找的分数总是存在的:
定理。给定包含至少一个分数的实数子区间J,J 包含唯一的最简单
分数。换句话说,有一个独特的分数f 使得
-
f 在 J 中,并且,
- 对于
J 中的任何其他分数g,f 比g 简单。
特别是,区间中最简单的分数总是有可能的最小分母,正如问题所要求的那样。定理中的“包含至少一个分数”条件是为了排除退化的情况,例如闭区间[√2, √2],它根本不包含任何分数。
我们的工作是编写一个函数,该函数接受有限浮点输入 x 并返回最简单的分数 n/d,其中 x 是最接近 n/d 的浮点数,采用目标浮点格式.
假设一个合理的浮点格式和舍入模式,四舍五入到x 的实数集将形成实线的非空子区间,具有合理的端点。所以我们的问题自然地分解为两个子问题:
第二个问题更有趣,对平台和语言细节的依赖更少;让我们先解决这个问题。
找出区间中最简单的分数
假设 IEEE 754 浮点格式和默认的舍入到偶数舍入模式,舍入到给定浮点数的间隔将是开放的或封闭的;对于其他舍入模式或格式,它可能是半开的(一端打开,另一端关闭)。所以在本节中,我们只关注开区间和闭区间,但适应半开区间并不难。
假设J 是具有有理端点的实线的非空子区间。为简单起见,我们假设J 是正实线的子区间。如果不是,那么它要么包含0——在这种情况下,0/1 是J 中的最简单分数——或者它是负实数线的子区间,我们可以求反,找到最简单的分数,然后取反。
那么下面给出了一个简单的递归算法,用于在J中找到最简单的分数:
- 如果J包含
1,那么1/1是J中的最简分数
- 否则,如果
J 是(0, 1) 的子区间,则J 中的最简分数为1/f,其中f 是1/J 中的最简分数。 (这直接来自“最简单”的定义。)
- 否则,
J 必须是(1, ∞) 的子区间,J 中最简单的分数是q + f,其中q 是最大整数,使得J - q 仍在正实数范围内,而f 是J - q 中最简单的分数。
最后一个陈述的证明草图:如果a / b 是J 中的最简分数,而c / d 是J - q 中的最简分数,则a / b 比或等于@987654379 @,并且c / d 简单于或等于(a - qb) / b。所以b <= d、a <= c + qd、d <= b和c <= a - qb,然后是b = d和a = c + qd,所以c / d = a / b - q。
在类似 Python 的伪代码中:
def simplest_in_interval(J):
# Simplest fraction in a subinterval J of the positive reals
if J < 1:
return 1 / simplest_in_interval(1/J)
else if 1 < J:
q = largest_integer_to_the_left_of(J)
return q + simplest_in_interval(J - q)
else:
# If we get here then J contains 1.
return 1
要查看算法必须始终终止并且不能进入无限循环,请注意每个反演步骤之后必须有一个J - q 步骤,并且每个J - q 步骤都会减少左右端点的分子的区间。具体来说,如果区间的端点是a/b和c/d,则和abs(a) + abs(c) + b + d是一个正整数,随着算法的进行而逐渐减小。
要将以上内容翻译成真正的 Python 代码,我们必须处理一些细节。首先,我们现在假设J 是一个闭区间;我们将适应下面的开放区间。
我们将通过端点left 和right 来表示我们的区间,这两个端点都是正的fraction.Fraction 实例。那么下面的Python代码实现了上述算法。
from fractions import Fraction
from math import ceil
def simplest_in_closed_interval(left, right):
"""
Simplest fraction in [left, right], assuming 0 < left <= right < ∞.
"""
if right < 1: # J ⊂ (0, 1)
return 1 / simplest_in_closed_interval(1 / right, 1 / left)
elif 1 < left: # J ⊂ (1, ∞):
q = ceil(left) - 1 # largest q satisfying q < left
return q + simplest_in_closed_interval(left - q, right - q)
else: # left <= 1 <= right, so 1 ∈ J
return Fraction(1)
这是一个运行示例:
>>> simplest_in_closed_interval(Fraction("3.14"), Fraction("3.15"))
Fraction(22, 7)
原则上,开区间的代码同样简单,但在实践中存在一个复杂问题:我们可能需要处理无限区间。例如,如果我们的原始区间是J = (2, 5/2),那么第一步将该区间移动2 得到(0, 1/2),然后将该区间倒转得到(2, ∞)。
因此,对于开区间,我们将继续通过其端点的一对 (left, right) 来表示我们的区间,但现在 right 要么是 fractions.Fraction 实例,要么是特殊常量 INFINITY。而不是简单地使用1 / left 来取左端点的倒数,我们需要一个辅助函数来计算分数或INFINITY 的倒数,以及另一个用于减法的辅助函数,确保INFINITY - q 提供INFINITY。以下是这些辅助功能:
#: Constant used to represent an unbounded interval
INFINITY = "infinity"
def reciprocal(f):
""" 1 / f, for f either a nonnegative fraction or ∞ """
if f == INFINITY:
return 0
elif f == 0:
return INFINITY
else:
return 1 / f
def shift(f, q):
""" f - q, for f either a nonnegative fraction or ∞ """
if f == INFINITY:
return INFINITY
else:
return f - q
这是主要功能。注意if 和elif 条件中不等式的变化,以及我们现在想使用floor(left) 而不是ceil(left) - 1 来找到位于区间左侧的最大整数q:
from fractions import Fraction
from math import floor
def simplest_in_open_interval(left, right):
"""
Simplest fraction in (left, right), assuming 0 <= left < right <= ∞.
"""
if 1 <= left: # J ⊆ (1, ∞):
q = floor(left)
return q + simplest_in_open_interval(shift(left, q), shift(right, q))
elif right != INFINITY and right <= 1: # J ⊆ (0, 1)
return 1 / simplest_in_open_interval(reciprocal(right), reciprocal(left))
else: # left < 1 < right, so 1 ∈ J
return Fraction(1)
上面的代码是为了清晰而不是效率而优化的:它在大复杂度方面相当有效,但在实现细节方面却不是。我把它留给读者转换成更有效的东西。第一步是使用整数分子和分母,而不是 fractions.Fraction 实例。如果您对它的外观感兴趣,请查看我在 PyPI 上的 simplefractions 包中的实现。
找到四舍五入到给定浮点数的区间
现在我们可以找到给定区间中最简单的分数,我们需要解决问题的另一半:找到四舍五入到给定浮点数的区间。执行此操作的细节在很大程度上取决于语言、使用的浮点格式,甚至是使用的舍入模式。
这里我们概述了一种在 Python 中执行此操作的方法,假设 IEEE 754 binary64 浮点格式具有默认的round-ties-to-even舍入模式。
为简单起见,假设我们的输入 float x 是正数(并且是有限的)。
Python >= 3.9 提供了一个函数math.nextafter,它允许我们从x 检索下一个上下浮动。下面是一个对最接近 π 的浮点数执行此操作的示例:
>>> import math
>>> x = 3.141592653589793
>>> x_plus = math.nextafter(x, math.inf)
>>> x_minus = math.nextafter(x, 0.0)
>>> x_plus, x_minus
(3.1415926535897936, 3.1415926535897927)
(请注意,通常要这样做,我们还需要处理x 是最大可表示浮点数并且math.nextafter(x, math.inf) 给出无穷大的特殊情况。)
四舍五入到x 的区间边界在x 和相邻浮点数之间。 Python 允许我们将浮点数转换为相应的精确值作为分数:
>>> from fractions import Fraction
>>> left = (Fraction(x) + Fraction(x_minus)) / 2
>>> right = (Fraction(x) + Fraction(x_plus)) / 2
>>> print(left, right)
14148475504056879/4503599627370496 14148475504056881/4503599627370496
我们还需要知道我们是否有一个封闭或开放的区间。我们可以查看位表示来弄清楚(这取决于浮点数的最低有效位是0 还是1),或者我们可以测试一下我们的区间端点是否舍入到x 或不是:
>>> float(left) == x
True
>>> float(right) == x
True
他们有,所以我们有一个封闭的区间。通过查看浮点数的十六进制表示可以证实这一点:
>>> x.hex()
'0x1.921fb54442d18p+1'
所以我们可以使用simplest_in_closed_interval 找到最简单的分数,该分数四舍五入为x:
>>> simplest_in_closed_interval(left, right)
Fraction(245850922, 78256779)
>>> 245850922 / 78256779 == x
True
把它们放在一起
虽然核心算法很简单,但有足够多的极端情况需要处理(负值、开区间与闭区间、sys.float_info.max 等),最终完整的解决方案有点过于混乱,无法完整发布这个答案。前段时间,我整理了一个名为simplefractions 的Python 包来处理所有这些极端情况;它是available on PyPI。它在行动中:
>>> from simplefractions import simplest_from_float
>>> simplest_from_float(0.1)
Fraction(1, 10)
>>> simplest_from_float(-3.3333333333333333)
Fraction(-10, 3)
>>> simplest_from_float(22/7)
Fraction(22, 7)
>>> import math
>>> simplest_from_float(math.pi)
Fraction(245850922, 78256779)