计算浮点值以将整数乘以产生小于 1 的最大数

Posted 2023-02-16

【中文标题】计算浮点值以将整数乘以产生小于 1 的最大数

【英文标题】：Calculating a floating-point value to multiply an integer with to produce largest number less than 1

【发布时间】：2019-07-30 14:41:54

【问题描述】：

对于给定的 32 位整数 y，选择一个双精度浮点数 x 以使 x*(double)y 尽可能大，但严格小于 1。

有没有算法可以计算出来？

我只是好奇，我没有无法通过其他方式解决的真正问题。

更新：

问题有点不准确，所以这里更新一下：

浮点格式为 IEEE-754 binary64 乘法的结果不是数学实数，而是浮点数

感谢@DavidEisenstat 和@EricPostpischil 注意到这些。

【问题讨论】：

浮点结果小于1或数学结果小于1？

双精度是指 IEEE-754 binary64 吗？

如果x 和y 的实数乘积为1−2^−54 或更大，则四舍五入为1（使用四舍五入时，绑定到-甚至）。因此，对于正数y，x 必须是小于 (1−2^−54)/y 的最大可表示数。

因为 y 被限制为 31 位（31 加一个符号，-2^31 单独求解 [琐碎]），如果 x•y 接近 1−2 ^−54，x 的 LSB 必须至少为 2^−84。所以，如果我们可以计算出商 (1−2^−54)/y 到那个精度，我们就可以截断后续位，这就是所需的x。

【参考方案1】：

这里有一个算法可以解决你的问题，但依赖于一些位旋转。

先计算并存储

x_attempt = 1.0 / (double) y

该值将非常接近您想要的 x 值，但由于浮点数学中的近似值，它可能会有些偏差。

现在检查x_attempt * (double) y 的值。如果该值不严格小于1.0，则将x_attempt 的值“微调”到下一个可以存储在双精度变量中的较小值。再次检查x_attempt * (double) y 并继续向下移动，直到该值严格小于1.0。这使用了一个循环，但它执行的次数很少，假设您的平台上的浮点运算完全没有问题。

如果该值严格小于1.0，则“微调”x_attempt 的值，直到该乘积为1.0 或更大。然后将x设置为x_attempt之前的值。

当那是我的首选语言时，我在 Borland Delphi 中实现了我自己的“轻推”例程。询问您在为您的语言和环境编写此类例程时是否需要帮助。

---

您的问题促使我用我当前的首选语言 Python 3.7 编写“微调”上下例程。它们在这里，没有我的单元测试例程。

import math

MIN_DOUBLE = 2 ** -1074  # 4.940656458412465442e-324
MIN_DOUBLE_DENORMAL = MIN_DOUBLE
MIN_DOUBLE_NORMAL = 2 ** -1022  # 2.225073858507201383e-308
MAX_DOUBLE = 1.7976931348623157082e308  # just under 2 ** 1024
EPSILON_DOUBLE_HALF = 2 ** -53  # changes most significand values
EPSILON_DOUBLE_FOURTH = 2 ** -54  # changes significand of -0.5


def nudged_up(x: float) -> float:
    """Return the next larger float value.

    NOTES:  1.  This routine expects that `float` values are implemented
                as  IEEE-754 binary64 and includes denormal values. No
                check is done on these assumptions.
            2.  This raises an OverflowError for MAX_DOUBLE. All other
                float values do not raise an error.
            3.  For +INF, -INF, or NAN, the same value is returned.
    """
    if x == 0.0:
        return MIN_DOUBLE_DENORMAL
    significand, exponent = math.frexp(x)
    if exponent < -1021:  # denormal
        return x + MIN_DOUBLE_DENORMAL
    if significand == -0.5:  # nudging will change order of magnitude
        diff = EPSILON_DOUBLE_FOURTH
    else:  # usual case
        diff = EPSILON_DOUBLE_HALF
    return math.ldexp(significand + diff, exponent)


def nudged_down(x: float) -> float:
    """Return the next smaller float value.

    NOTES:  1.  This routine expects that `float` values are implemented
                as  IEEE-754 binary64 and includes denormal values. No
                check is done on these assumptions.
            2.  This raises an OverflowError for -MAX_DOUBLE. All other
                float values do not raise an error.
            3.  For +INF, -INF, or NAN, the same value is returned.
    """
    return -nudged_up(-x)

这里是 Python 3.7 代码，可以回答您的问题。请注意，如果输入参数为零，则会引发错误——您可能需要更改它。

def lower_reciprocal(y: int) -> float:
    """Given a (32-bit) integer y, return a float x for which
    x * float(y) is as large as possible but strictly less than 1.

    NOTES:  1.  If y is zero, a ZeroDivisionError exception is raised.
    """
    if y < 0:
        return -lower_reciprocal(-y)
    float_y = float(y)
    x = 1.0 / float_y
    while x * float_y < 1.0:
        x = nudged_up(x)
    while x * float_y >= 1.0:
        x = nudged_down(x)
    return x

【参考方案2】：

这可能不是您的想法，但回答此类问题的一个很好的选择是使用现成的 SMT 求解器，它可以解决许多领域的约束满足问题；包括 IEEE 浮点数。详情请见http://smtlib.cs.uiowa.edu/。

Microsoft 的 Z3 是一款优秀（而且免费！）的 SMT 求解器：https://github.com/Z3Prover/z3

SMT 求解器使用所谓的 SMTLib 语言编写脚本，该语言更面向机器。但它们也为许多语言（C、C++、Java、Python）提供编程 API，并且有许多高级语言绑定使它们易于使用，包括来自 Scala、O'Caml、Go、Haskell 的那些很少。

例如，下面是如何在 Haskell 的 SBV 绑定到 Z3 中对查询进行编码。如果您不阅读 Haskell，请不要担心；这个想法是为了说明如何在非常高的水平上快速编写此类问题：

import Data.SBV

q :: SInt32 -> IO Double
q y = do LexicographicResult m <- optimize Lexicographic $ do
                                      x <- sDouble "x"
                                      let res = (sFromIntegral (y::SInt32) * x)
                                      constrain $ res .< 1
                                      maximize "mx" res
         case getModelValue "x" m of
           Just x  -> return x
           Nothing -> error "No such value exists!" -- shouldn't happen!

通过他的一段代码，我们可以使用 Haskell 解释器ghci (https://www.haskell.org/ghc/) 来查询各种值：

*Main> q 12
8.333333333333331e-2
*Main> q 821
1.2180267965895247e-3

这些答案可以快速计算，并且可以进一步使用并与其他约束相结合，以轻松找到此类感兴趣的值。

如果您有兴趣，还可以查看 z3 的 Python 绑定，这可能是最容易上手的，尽管不像其他绑定那样类型安全，正如您想象的那样。见这里：https://ericpony.github.io/z3py-tutorial/guide-examples.htm