Python乘法失去浮点精度

Posted 2023-02-16

【中文标题】Python乘法失去浮点精度

【英文标题】：Python multiplication losing floating point precision

【发布时间】：2020-03-01 04:47:59

【问题描述】：

我最近在处理浮点数，我意识到浮点数有一些我没想到的东西。这是一个例子

a = 0.1
print(f"a:0.20f")
#'0.10000000000000000555'
b = a * 10
print(f"b:0.20f")
#'1.00000000000000000000'

我希望最后一次打印输出1.00000000000000005551（即，1 后跟0.1 的数字 1 到 21）。

我很好奇的是为什么浮点误差在乘以 10 时会消失。正常的算术规则表明浮点误差会被传播，但实际上并没有发生。为什么会发生这种情况？有没有办法避免它？

【问题讨论】：

这能回答你的问题吗？ Is floating point math broken?

没有。我显然知道小数可能有浮点错误。我很好奇的是为什么浮点误差在乘以 10 时会消失。正常的算术规则表明浮点误差会被传播，但实际上并没有发生。

【参考方案1】：

10和0.1000000000000000055511151231257827021181583404541015625，该IEEE 754的64位的0.1二进制表示确切的实数算术乘积，是1.000000000000000055511151231257827021181583404541015625。 P>

它不能完全表示。它被 1.0 和 1.0000000000000002220446049250313080847263336181640625 括起来

它更接近 1.0，所以这是乘法的四舍五入结果。

我使用一个简短的 Java 程序计算了这些数字：

import java.math.BigDecimal;

public strictfp class Test 
    public static void main(String[] args) 
        BigDecimal rawTenth = new BigDecimal(0.1);
        BigDecimal realProduct = rawTenth.multiply(BigDecimal.TEN);
        System.out.println(realProduct);
        System.out.println(new BigDecimal(Math.nextUp(1.0)));
    

输出：

1.0000000000000000555111512312578270211815834045410156250
1.0000000000000002220446049250313080847263336181640625

【讨论】：

所以浮点数有一些界限，低于这个界限 Python 会在乘法后自动向下取整？是否有相关文档或标准？

不，不管是乘法、除法、加法等，还是数字的大小。它在round half even 模式下运行，并进行了计算，其实数结果比任何其他可表示数字更接近 1.0。

相关标准是 IEEE 754。该标准本身在付费墙后面，但有一篇有用的***文章 IEEE 754。

另一种消除乘法特殊处理理论的方法是使用 1.000000000000000055511151231257827021181583404541015625 作为文字并打印出来。

【参考方案2】：

这个答案展示了如何确定将 1/10 转换为浮点数并乘以 10 将只使用一点算术运算就可以准确地产生 1；无需计算大的或精确的数字。

您的 Python 实现使用常见的 IEEE-754 binary64 格式。 （Python 对应该使用哪种浮点格式实现并不严格。）在这种格式中，数字实际上表示为应用于某个 53 位整数乘以 2 的某个幂的符号（+ 或 -）。因为 2⁻⁴ ≤ 1/10 −3，所以最接近 1/10 的可表示数是某个整数 M 乘以 2^{−3−53支持>。 （-53 将 53 位整数缩放到 ½ 和 1 之间，-3 将其缩放到 2-4 和 2-3 之间。）我们称之为可表示的数 x。}

那么我们有 x = M•2⁻⁵⁶ = 1/10 + e，其中 e 是当我们将 1/10 舍入到最接近的可表示值时出现的一些舍入误差。由于我们四舍五入到最接近的可表示值，|e| ≤ ½•2^-56 = 2^-57.

要准确找出 e 是什么，请将 1/10 乘以 2⁵⁶。 WolframAlpha 告诉我们它是 7205759403792793+3/5。为了得到最接近的可表示值，我们应该四舍五入，所以 M = 7205759403792794 和 e = 2/5 • 2⁻⁵⁶。虽然我用 WolframAlpha 来说明这一点，但我们不需要 M，我们可以通过观察两个模 10 的幂的模式来找到 e：2¹→2, 2²→4, 2³→8, 2⁴→6, 2⁵→2, 2⁶→ 4，因此模式以 4 为循环重复，56 模 4 为 0，所以 2⁵⁶ 模 10 与 2⁴ 具有相同的余数，6，所以分数是 6/10 = 3/5。我们知道应该四舍五入到最接近的整数 1，所以 e = 2/5 • 2⁻⁵⁶。

所以 x = M•2^-56 = 1/10 + 2/5•2^-56。

现在我们可以算出用浮点运算计算 10•x 的结果了。结果就像我们首先用实数算术计算 10•x，然后四舍五入到最接近的可表示值。在实数算术中，10•x = 10•(1/10 + 2/5•2⁻⁵⁶) = 1 + 10•2/5•2⁻⁵⁶ = 1 + 4•2^-56 = 1 + 2^-54。两个相邻的可表示值是 1 和 1 + 2^-52，1 + 2^-54 比 1 + 2^{-52 更接近 1} 。所以结果是 1。

【讨论】：

我在注释格式方面遇到了问题，但我得到了 1.0 和实数乘积之间的差异作为 Math.pow(2,-54)，而不是 Math.pow(2,-53)。

@PatriciaShanahan：是的，我将 4•2^-56 错误地视为 2^-53 而不是 2^-54。已修复，谢谢。

为了比较，我在答案中添加了我用来计算数字的方法。尽管作为一名数学家受过早期培训，但我发现一点 Java hacking 更容易。