将浮点数转换为保证转换回原始浮点数的有理数
Convert a float to a rational number that is guaranteed to convert back to the original float
我正在寻找一种将浮点数转换为有理数的算法，这样可以保证有理数的计算结果返回到原始浮点数，并且分母最小。
一个朴素的算法可以直接return浮点数的实际值作为X/2N，但是2N 对于不是有限二进制分数的任何东西来说往往都非常大。例如，当数字0.1存储在双精度浮子中时，实际上是由³⁶⁰²⁸⁷⁹⁷⁰²⁸⁷⁹⁷⁰⁄₃₆₀₂₈₇₉₇₀₁₈₉₆₃₉₆₈（分母为2 55 ）的近似。但是，将 0.1 转换为 ¹⁄₁₀ 显然更好，而 ¹⁄₁₀ 将计算为 ³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹7⁄₃₆₀₈₇₉₇₀₁₈₉=2 under arithmetic=1₆₈₉=1
一个相关的问题是用最少的数字打印十进制浮点数（这个paper描述了一些技术），这可以被认为是这个问题的一个特殊版本，有一个额外的限制，即分母必须是10 的幂。
有 an existing questions，可能更多，但它们没有转换后的有理数必须计算为原始浮点数的约束。
让我们从一个定义开始，该定义准确地确定了在任何特定情况下我们正在寻找的分数：
Definition. Say that a fraction a/b is simpler than another fraction c/d (with both fractions written in lowest terms, with
positive denominators) if b <= d, abs(a) <= abs(c), and at least
one of those two inequalities is strict.
所以例如5/7比6/7简单，5/7比5/8简单，但是2/5和3/4都不是比另一个简单。 （我们这里没有总排序。）
那么根据这个定义，有一个不是很明显的定理，它保证我们要找的分数总是存在的：
Theorem. Given a subinterval J of the real numbers that contains at least one fraction, J contains a unique simplest
fraction. In other words, there's a unique fraction f such that
f is in J and,- for any other fraction
g in J, f is simpler than g.
特别是，区间中最简单的分数将始终具有尽可能小的分母，如问题中所要求的那样。定理中的“至少包含一个分数”条件是为了排除退化情况，例如闭区间 [√2, √2]，它根本不包含任何分数。
我们的工作是编写一个函数，该函数采用有限浮点输入 x 和 returns 最简单的分数 n/d 其中 x 是最接近的浮点数到 n/d，目标浮点格式。
假设合理合理的浮点格式和舍入模式，舍入到 x 的实数集将形成实线的非空子区间，具有有理端点。所以我们的问题自然分解为两个子问题：
问题1.给定一个目标浮点格式的浮点数x，描述所有舍入到[=的值的区间41=] 根据该浮点格式的规则。这涉及识别该区间的端点并确定该区间是开放的、封闭的还是半开放的。
问题 2.给定实线的有理端点的非空子区间 J，计算该子区间中的最简单分数。
第二个问题更有趣，对平台和语言细节的依赖更少；让我们先解决这个问题。
寻找区间中最简单的分数
假定 IEEE 754 浮点格式和默认的舍入到偶数舍入模式，给定浮点数的区间舍入将是开放式或封闭式；对于其他舍入模式或格式，它可能是半开的（一端打开，另一端关闭）。所以本节我们只看开区间和闭区间，但适应半开区间并不难。
假设J是实线的有理端点的非空子区间。为简单起见，我们假设 J 是 positive 实数线的子区间。如果不是，那么它要么包含 0 — 在这种情况下 0/1 是 J 中最简单的分数 — 或者它是负实数线的子区间，我们可以取反，找到最简单的分数，然后取反。
然后下面给出了一个简单的递归算法，用于在J中找到最简单的分数：
- 如果J包含
1，则1/1是J 中最简单的分数
- 否则，如果
J是(0, 1)的子区间，那么J中的最简分数是1/f，其中f是最简分数在 1/J。 （这直接来自 'simplest' 的定义。）
- 否则，
J一定是(1, ∞)的子区间，J中的最简分数是q + f，其中q是最大整数这样 J - q 仍然在正实数范围内，并且 f 是 J - q. 中最简单的分数
最后一个陈述的证明草图：如果 a / b 是 J 中最简单的分数并且 c / d 是 J - q 中最简单的分数，那么 a / b 比 (c + qd) / d 简单或等于，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)
上面的代码是为了清晰而不是效率而优化的：它在大 Oh 复杂性方面相当高效，但在实现细节方面却不是。我把它留给 reader 来转换成更有效的东西。第一步是始终使用整数分子和分母，而不是 fractions.Fraction 个实例。如果您对它的外观感兴趣，请查看我在 PyPI 上的 simplefractions 包中的实现。
查找舍入到给定浮点数的区间
现在我们可以找到给定区间内的最简单分数，我们需要解决问题的另一半：找到给定浮点数四舍五入的区间。这样做的细节将在很大程度上取决于语言、使用的浮点格式，甚至是使用的舍入模式等因素。
这里我们概述了在 Python 中执行此操作的一种方法，假设 IEEE 754 二进制 64 浮点格式具有默认的舍入到偶数舍入模式。
为简单起见，假设我们的输入浮点数 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 等），以至于一个完整的解决方案最终会有点过于混乱而无法解决post 在此答案中完整。前段时间我整理了一个名为 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)
我正在寻找一种将浮点数转换为有理数的算法，这样可以保证有理数的计算结果返回到原始浮点数，并且分母最小。
一个朴素的算法可以直接return浮点数的实际值作为X/2N，但是2N 对于不是有限二进制分数的任何东西来说往往都非常大。例如，当数字0.1存储在双精度浮子中时，实际上是由³⁶⁰²⁸⁷⁹⁷⁰²⁸⁷⁹⁷⁰⁄₃₆₀₂₈₇₉₇₀₁₈₉₆₃₉₆₈（分母为2 55 ）的近似。但是，将 0.1 转换为 ¹⁄₁₀ 显然更好，而 ¹⁄₁₀ 将计算为 ³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹7⁄₃₆₀₈₇₉₇₀₁₈₉=2 under arithmetic=1₆₈₉=1
一个相关的问题是用最少的数字打印十进制浮点数（这个paper描述了一些技术），这可以被认为是这个问题的一个特殊版本，有一个额外的限制，即分母必须是10 的幂。
有 an existing questions，可能更多，但它们没有转换后的有理数必须计算为原始浮点数的约束。
让我们从一个定义开始，该定义准确地确定了在任何特定情况下我们正在寻找的分数：
Definition. Say that a fraction
a/bis simpler than another fractionc/d(with both fractions written in lowest terms, with positive denominators) ifb <= d,abs(a) <= abs(c), and at least one of those two inequalities is strict.
所以例如5/7比6/7简单，5/7比5/8简单，但是2/5和3/4都不是比另一个简单。 （我们这里没有总排序。）
那么根据这个定义，有一个不是很明显的定理，它保证我们要找的分数总是存在的：
Theorem. Given a subinterval
Jof the real numbers that contains at least one fraction,Jcontains a unique simplest fraction. In other words, there's a unique fractionfsuch that
fis inJand,- for any other fraction
ginJ,fis simpler thang.
特别是，区间中最简单的分数将始终具有尽可能小的分母，如问题中所要求的那样。定理中的“至少包含一个分数”条件是为了排除退化情况，例如闭区间 [√2, √2]，它根本不包含任何分数。
我们的工作是编写一个函数，该函数采用有限浮点输入 x 和 returns 最简单的分数 n/d 其中 x 是最接近的浮点数到 n/d，目标浮点格式。
假设合理合理的浮点格式和舍入模式，舍入到 x 的实数集将形成实线的非空子区间，具有有理端点。所以我们的问题自然分解为两个子问题：
问题1.给定一个目标浮点格式的浮点数
x，描述所有舍入到[=的值的区间41=] 根据该浮点格式的规则。这涉及识别该区间的端点并确定该区间是开放的、封闭的还是半开放的。问题 2.给定实线的有理端点的非空子区间
J，计算该子区间中的最简单分数。
第二个问题更有趣，对平台和语言细节的依赖更少；让我们先解决这个问题。
寻找区间中最简单的分数
假定 IEEE 754 浮点格式和默认的舍入到偶数舍入模式，给定浮点数的区间舍入将是开放式或封闭式；对于其他舍入模式或格式，它可能是半开的（一端打开，另一端关闭）。所以本节我们只看开区间和闭区间，但适应半开区间并不难。
假设J是实线的有理端点的非空子区间。为简单起见，我们假设 J 是 positive 实数线的子区间。如果不是，那么它要么包含 0 — 在这种情况下 0/1 是 J 中最简单的分数 — 或者它是负实数线的子区间，我们可以取反，找到最简单的分数，然后取反。
然后下面给出了一个简单的递归算法，用于在J中找到最简单的分数：
- 如果J包含
1，则1/1是J中最简单的分数
- 否则，如果
J是(0, 1)的子区间，那么J中的最简分数是1/f，其中f是最简分数在1/J。 （这直接来自 'simplest' 的定义。） - 否则，
J一定是(1, ∞)的子区间，J中的最简分数是q + f，其中q是最大整数这样J - q仍然在正实数范围内，并且f是J - q. 中最简单的分数
最后一个陈述的证明草图：如果 a / b 是 J 中最简单的分数并且 c / d 是 J - q 中最简单的分数，那么 a / b 比 (c + qd) / d 简单或等于，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)
上面的代码是为了清晰而不是效率而优化的：它在大 Oh 复杂性方面相当高效，但在实现细节方面却不是。我把它留给 reader 来转换成更有效的东西。第一步是始终使用整数分子和分母，而不是 fractions.Fraction 个实例。如果您对它的外观感兴趣，请查看我在 PyPI 上的 simplefractions 包中的实现。
查找舍入到给定浮点数的区间
现在我们可以找到给定区间内的最简单分数，我们需要解决问题的另一半：找到给定浮点数四舍五入的区间。这样做的细节将在很大程度上取决于语言、使用的浮点格式，甚至是使用的舍入模式等因素。
这里我们概述了在 Python 中执行此操作的一种方法，假设 IEEE 754 二进制 64 浮点格式具有默认的舍入到偶数舍入模式。
为简单起见，假设我们的输入浮点数 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 等），以至于一个完整的解决方案最终会有点过于混乱而无法解决post 在此答案中完整。前段时间我整理了一个名为 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)