问题描述
我正在寻找一种将浮点数转换为有理数的算法,以便保证有理数计算回原始浮点数,并最小化分母。
一个简单的算法只能将浮点数的实际值返回为 X / 2N,但是对于任何不是有限二进制分数。例如,当数字 0.1 存储在双精度浮点数中时,实际上近似为但是,将 0.1 转换为 ¹⁄₁₀ 显然更好,并且¹⁄₁₀ 将计算为 ³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹⁷⁄₃₆₀₂₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₆₇₆浮点
一个相关的问题是用最少的数字打印十进制浮点数(这个 paper 描述了一些技术),这可以被认为是这个问题的一个特殊版本,有一个额外的约束,即分母必须是一个幂共 10 个。
有 an existing questions,而且可能更多,但它们没有转换后的有理数必须计算为原始浮点数的约束。
解决方法
让我们从一个定义开始,该定义准确地确定我们在任何特定情况下寻找的分数:
定义。 假设一个分数 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
中,并且, - 对于
g
中的任何其他分数J
,f
比g
简单。
特别是,根据问题的要求,区间中最简单的分数总是具有尽可能小的分母。定理中的“至少包含一个分数”条件是为了排除像闭区间 [√2,√2]
这样根本不包含任何分数的退化情况。
我们的工作是编写一个函数,该函数接受有限浮点输入 x
并返回最简单的分数 n/d
,其中 x
是最接近 n/d
的浮点数,目标浮点格式。
假设合理合理的浮点格式和舍入模式,舍入到 x
的实数集将形成实数线的非空子区间,具有合理的端点。所以我们的问题自然分解为两个子问题:
-
问题 1. 给定目标浮点格式的浮点数
x
,根据以下规则描述四舍五入为x
的所有值的区间那种浮点格式。这涉及识别该区间的端点并确定该区间是开放、封闭还是半开放。 -
问题 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}}。 (这是从“最简单”的定义中直接得出的。) - 否则,
1/J
必须是J
的子区间,(1,∞)
中最简单的分数是J
,其中q + f
是最大的整数,使得q
仍在正实数内,J - q
是f
中最简单的分数。
对于最后一个陈述的证明草图:如果 J - q
是 a / b
中的最简单分数,而 J
是 c / d
中的最简单分数,那么 {{1 }} 比J - q
更简单,a / b
比(c + qd) / d
更简单。所以c / d
、(a - qb) / b
、b <= d
和a <= c + qd
,然后是d <= b
和c <= a - qb
,所以b = d
。>
在类似 Python 的伪代码中:
a = c + qd
要看到算法必须总是终止并且不能进入无限循环,注意每一个求逆步骤后面必须跟一个c / d = a / b - q
步骤,并且每一个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
是一个正整数,随着算法的进行而稳步减少。
要将上述内容翻译成真正的 Python 代码,我们必须处理一些细节。首先,让我们暂时假设 c/d
是一个闭区间;我们将适应下面的开放间隔。
我们将用端点 abs(a) + abs(c) + b + d
和 J
表示我们的区间,这两个端点都是正 left
实例。然后下面的Python代码实现了上面的算法。
right
这是一个运行示例:
fraction.Fraction
原则上,开区间的代码同样简单,但实际上有一个复杂之处:我们可能需要处理无限区间。例如,如果我们的原始区间是 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)
实例,要么是一个特殊的常量 {{ 1}}。而不是简单地能够使用 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
来找到最大整数 {{1 }} 位于区间左侧:
elif
上面的代码是为了清晰而不是效率而优化的:它在 big-Oh 复杂性方面相当有效,但在实现细节方面却没有。我把它留给读者来转换成更有效的东西。第一步是使用整数分子和分母而不是 floor(left)
实例。如果您对它的外观感兴趣,请查看我在 PyPI 上的 simplefractions 包中的实现。
找到舍入到给定浮点数的间隔
既然我们可以在给定的区间内找到最简单的分数,我们需要解决问题的另一半:找到四舍五入到给定浮点数的区间。这样做的细节将在很大程度上取决于语言、使用的浮点格式,甚至使用的舍入模式等。
这里我们概述了在 Python 中执行此操作的一种方法,假设 IEEE 754 binary64 浮点格式具有默认的舍入到偶数舍入模式。
为简单起见,假设我们的输入浮点数 ceil(left) - 1
是正的(并且是有限的)。
Python >= 3.9 提供了一个函数 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
(注意,一般来说,要做到这一点,我们还需要处理特殊情况,其中 x
是最大的可表示浮点数,而 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)
和相邻浮点数之间。 Python 允许我们将浮点数转换为相应的精确值作为分数:
x
我们还需要知道我们有一个闭区间还是一个开区间。我们可以查看位表示来弄清楚(这取决于浮点数的最低有效位是 math.nextafter(x,math.inf)
还是 x
),或者我们可以测试一下我们的区间端点是否四舍五入到x
与否:
>>> 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'
等),以至于一个完整的解决方案最终有点过于混乱而无法发布全在这个答案。前段时间我整理了一个名为 x
的 Python 包来处理所有这些极端情况;它是available on PyPI。这是在行动:
simplest_in_closed_interval