15

Python 是否提供了一个函数来获取通过增加现有浮点值的最低有效位而产生的浮点值?

我正在寻找类似于std::nextafterC++11 中添加的函数的东西。

4

3 回答 3

26

这里有五个(实际上是四个半)可能的解决方案。

解决方案一:使用 Python 3.9 或更高版本

2020 年 10 月发布的 Python 3.9 包含一个新的标准库函数math.nextafter,该函数直接提供此功能:用于math.nextafter(x, math.inf)获取下一个正无穷大的浮点数。例如:

>>> from math import nextafter, inf
>>> nextafter(100.0, inf)
100.00000000000001

如果您查看该方法提供的十六进制表示,则更容易验证此函数是否确实产生float.hex了下一个浮动:

>>> 100.0.hex()
'0x1.9000000000000p+6'
>>> nextafter(100.0, inf).hex()
'0x1.9000000000001p+6'

Python 3.9 还引入了一个密切相关且经常有用的伴随函数math.ulp,它给出了一个值和下一个远离零的值之间的差异:

>>> from math import ulp
>>> nextafter(100.0, inf) - 100.0
1.4210854715202004e-14
>>> ulp(100.0)
1.4210854715202004e-14

解决方案 2:使用 NumPy

如果您没有 Python 3.9 或更高版本,但您可以访问 NumPy,那么您可以使用numpy.nextafter. 对于常规float的 Python,语义匹配math.nextafter(尽管说 Python 的语义匹配 NumPy 更公平,因为 NumPy 早在 Python 之前就提供了这个功能

>>> from numpy import nextafter, inf
>>> nextafter(100.0, inf)
100.00000000000001

nextafter解决方案3:自己包装C

C 指定一个nextafter函数math.h(例如参见 C99 的第 7.12.11.3 节);这正是 Python >= 3.9 在其math模块中封装和公开的函数。如果您没有 Python 3.9 或更高版本,您可以使用ctypescffi动态调用 C's nextafter,或者编写一个简单的Cython包装器或 Python C 扩展来公开 C's nextafter。如何做到这一点的细节已经在其他地方得到了很好的解释:在@ Endophage对这个问题的回答中,以及在这个对类似 StackOverflow 问题的回答中(这个问题作为副本被关闭)。

解决方案 4:通过struct模块进行位操作

如果您愿意(在实践中几乎总是安全的)假设 Python 使用 IEEE 754 浮点,那么编写一个 Python 函数来提供nextafter. 需要一点点小心才能正确处理所有角落案例。

IEEE 754 二进制浮点格式经过巧妙设计,因此从一个浮点数移动到“下一个”浮点数就像增加位表示一样简单。这适用于 range 中的任何数字[0, infinity),正好跨越指数边界和次正规。要生成一个nextUp涵盖整个浮点范围的版本,您还需要处理负数、无穷大、nans 和一个涉及负零的特殊情况。nextUp以下是Python中 IEEE 754 函数的标准兼容版本。它涵盖了所有极端情况。

import math
import struct

def nextup(x):
    # NaNs and positive infinity map to themselves.
    if math.isnan(x) or (math.isinf(x) and x > 0):
        return x

    # 0.0 and -0.0 both map to the smallest +ve float.
    if x == 0.0:
        x = 0.0

    n = struct.unpack('<q', struct.pack('<d', x))[0]
    if n >= 0:
        n += 1
    else:
        n -= 1
    return struct.unpack('<d', struct.pack('<q', n))[0]

nextDown然后的实现nextAfter看起来像这样。(请注意,这nextAfter不是 IEEE 754 指定的函数,因此对于 IEEE 特殊值应该发生什么有一些猜测。这里我遵循 Pythondecimal.Decimal类所基于的 IBM 十进制算术标准。)

def nextdown(x):
    return -nextup(-x)

def nextafter(x, y):
    # If either argument is a NaN, return that argument.
    # This matches the implementation in decimal.Decimal
    if math.isnan(x):
        return x
    if math.isnan(y):
        return y

    if y == x:
        return y
    elif y > x:
        return nextup(x)
    else:
        return nextdown(x)

(部分)解决方案5:浮点运算

如果x是一个积极的不太小的float,并且你愿意假设 IEEE 754 binary64 格式和语义,有一个非常简单的解决方案:下一个从xis向上浮动,下一个从isx / (1 - 2**-53)向下浮动。xx * (1 - 2**-53)

更详细地说,假设以下所有条件都是正确的:

  • 您不关心 IEEE 754 极端情况(零、无穷大、次正规、nans)
  • 您不仅可以假设 IEEE 754 binary64 浮点格式,还可以假设 IEEE 754 binary64语义:即所有基本算术运算都根据当前舍入模式正确舍入
  • 您可以进一步假设当前舍入模式是 IEEE 754 默认的round-ties-to-even 模式。

那么这个数量1 - 2**-53可以精确地表示为 a float,并且给定一个正的非正常 Python float xx / (1 - 2**-53)将匹配nextafter(x, inf)。类似地,x * (1 - 2**-53)将匹配,除了在最小正正常值nextafter(x, -inf)的极端情况下, 。x2**-1022

使用它时要注意一件事:表达式2**-53将从系统的数学库中调用您,并且期望正确舍入pow通常是不安全的。pow有许多更安全的方法来计算这个常数,其中之一是使用float.fromhex. 这是一个例子:

>>> d = float.fromhex('0x1.fffffffffffffp-1')  # 1 - 2**-53, safely
>>> d
0.9999999999999999
>>> x = 100.0
>>> x / d  # nextup(x), or nextafter(x, inf)
100.00000000000001
>>> x * d  # nextdown(x), or nextafter(x, -inf)
99.99999999999999

这些技巧适用于正常的浮点数范围,包括像 2 的精确幂这样的尴尬情况。

对于一个证明的草图:为了显示与正法x / d线相匹配,我们可以按 2 的幂进行缩放而不影响正确性,因此在证明中我们可以假设 不失一般性。如果我们写出(被认为是实数,而不是浮点数)的精确数学值,则等于。结合不等式,我们可以得出结论,由于浮点数在区间内完全分开,足以保证最接近的浮点数是。的证明类似。nextafter(x, inf)x0.5 <= x < 1.0zx / dz - xx * 2**-53 / (1 - 2**-53)0.5 <= x <= 1 - 2**-532**-54 < z - x <= 2**-532**-53[0.5, 1.0]znextafter(x, inf)x * d

于 2012-05-03T06:10:52.267 回答
2

更新:

原来这是一个重复的问题(在 google 中作为搜索“c++ nextafter python”的结果 #2 出现):Increment a python floating point value by the minimum possible amount

公认的答案提供了一些可靠的解决方案。

原始答案:

当然这不是完美的解决方案,但使用 cython 只需几行就可以让您包装现有的 C++ 函数并在 Python 中使用它。我已经编译了下面的代码,它可以在我的 ubuntu 11.10 机器上运行。

首先,一个 .pyx 文件(我称为我的 nextafter.pyx)定义了您与 C++ 的接口:

cdef extern from "cmath":
    float nextafter(float start, float to)

def pynextafter(start, to):
    cdef float float_start = float(start)
    cdef float float_to = float(to)
    result = nextafter(start, to)
    return result

然后 setup.py 定义了如何构建扩展:

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext 

ext_modules=[
    Extension("nextafter",
        ["nextafter.pyx"],
        libraries=[],
        library_dirs=[],
        include_dirs=[],
        language="c++",
    )
]

setup(
    name = "nextafter",
    cmdclass = {"build_ext": build_ext},
    ext_modules = ext_modules
)

确保它们在同一目录中,然后使用python setup.py build_ext --inplace. 我希望你能看到如何将 nextafter 的其他变体添加到扩展中(对于双打等)。一旦建成,你应该有一个 nextafter.so。在同一目录中启动 python(或将 nextafter.so 放在您的路径上),您应该可以调用from nextafter import pynextafter.

享受!

于 2012-05-02T20:19:23.543 回答
0

查看http://docs.python.org/library/stdtypes.html#float.hex

让我们试试这个不太了解 next after 的实现。

首先,我们需要从十六进制字符串中提取十六进制部分和指数:

def extract_parts(hex_val):
    if not hex_val.startswith('0x1.'):
        return None
    relevant_chars = hex_val[4:]
    if not len(relevant_chars) > 14 and relevant_chars[13] == 'p':
        return None
    hex_portion = int(relevant_chars[:13], 16)
    if relevant_chars[14] == '+':
        p_val = int(relevant_chars[15:])
    elif relevant_chars[14] == '-':
        p_val = -int(relevant_chars[15:])
    else:
        return None
    return (hex_portion, p_val)

然后我们需要一种在正方向或负方向上递增的方法(我们假设十六进制字符串已经转换为整数hex_portion):

def increment_hex(hex_portion, p_val, direction):
    if hex_portion == 0 and direction == -1:
        new_hex = 'f' * 13
        p_val -= 1
    elif hex_portion == int('f' * 13, 16) and direction == 1:
        new_hex = '0' * 13
        p_val += 1
    else:
        new_hex = hex(hex_portion + direction)[2:].rstrip('L').zfill(13)

    if len(new_hex) != 13:
        return None
    return format_hex(new_hex, p_val)

我们需要一个辅助函数来重新格式化可接受的十六​​进制字符串和指数,我在上面使用过:

def format_hex(hex_as_str, p_val):
    sign = '-' if p_val < 0 else '+'
    return '0x1.%sp%s%d' % (hex_as_str, sign, p_val)

最后,实施nextafter

def nextafter(float_val):
    hex_equivalent = float_val.hex()
    hex_portion, p_val = extract_parts(hex_equivalent)

    direction = 1
    new_hex_equiv = increment_hex(hex_portion, p_val, direction)
    return float.fromhex(new_hex_equiv)
于 2012-05-02T20:09:22.213 回答