logo标签列表

浮点数和 decimal.Decimal 的小数位问题

pythonfloating-pointdecimalfloating-accuracy
3

我似乎在使用浮点数时失去了很多精度。

例如，我需要解决一个矩阵：

4.0x -2.0y 1.0z =11.0
1.0x +5.0y -3.0z =-6.0
2.0x +2.0y +5.0z =7.0

这是我用来从文本文件导入矩阵的代码：
f = open('gauss.dat')
lines =  f.readlines()
f.close()

j=0
for line in lines:
    bits = string.split(line, ',')
    s=[]
    for i in range(len(bits)):
        if (i!= len(bits)-1):
            s.append(float(bits[i]))
            #print s[i]
    b.append(s)
    y.append(float(bits[len(bits)-1]))

我需要使用高斯-塞德尔法求解，因此我需要重新排列x、y和z的方程：

x=(11+2y-1z)/4
y=(-6-x+3z)/5
z=(7-2x-2y)/7

我使用以下代码对方程进行重新排列。 b是系数矩阵，y是答案向量:

def equations(b,y):
    i=0
    eqn=[]
    row=[]
    while(i<len(b)):
        j=0
        row=[]
        while(j<len(b)):
            if(i==j):
                row.append(y[i]/b[i][i])
            else:
                row.append(-b[i][j]/b[i][i])
            j=j+1
        eqn.append(row)
        i=i+1
    return eqn

然而，我得到的答案并非精确到小数点。

例如，根据上述第二个方程式重新排列后，我应该得到：

y=-1.2-.2x+.6z

我得到的是：
y=-1.2-0.20000000000000001x+0.59999999999999998z

这可能看起来不是一个大问题，但是当你把数字提高到非常高的幂时，误差就会相当大。有什么解决方法吗？我尝试了使用Decimal类，但它在幂运算方面表现不佳（例如，Decimal(x)**2）。 有什么想法吗？
6个回答
14

IEEE浮点数是二进制的，而不是十进制的。没有固定长度的二进制小数可以准确地表示0.1或任何倍数。它是一个循环小数，就像十进制中的1/3。

请阅读《计算机科学家应该了解的浮点算术知识》

除了Decimal类之外，其他选择有：

  • 使用Common Lisp或Python 2.6或另一种精确分数运算的语言

  • 将double转换为接近的分数，例如使用frap

@Doug：我添加了对Python 2.6及其fractions模块的引用。 - tzot
只需使用类似gmpy的有理数算术库即可。 - Brian
1
你是什么意思说“任何倍数”？有一个固定长度的二进制小数，恰好为0.5。 - Nosredna
是的，但是没有固定长度的二进制小数可以精确表示0.1（十分之一）或任何0.1的倍数。 - Doug Currie
@Doug Currie 0.5是0.1的倍数，所以最后一部分不正确 ;) - Voo
12

我对 Decimal 类不够熟悉，无法帮助你解决问题，但你的问题是由于十进制小数在二进制中往往无法以准确方式表示，因此你看到的是最接近的可能近似值。除非使用特殊类（例如 Decimal，可能性比较大），否则无法避免这个问题。

编辑：关于 Decimal 类，什么没有正常工作？只要从字符串开始，而不是从浮点数开始，幂运算似乎都能正常工作。

>>> import decimal
>>> print(decimal.Decimal("1.2") ** 2)
1.44

模块文档相当清晰地解释了为什么需要使用decimal.Decimal以及其使用方法，如果您还没有查看过，请务必查看。

我之前遇到的关于幂的问题是，我试图将其提高为0.5次方。为此，我必须编写decimal.Decimal("1.2") ** decimal.Decimal("0.5")。 - darudude
你缺少一个引号。也许那就是问题所在。如果你觉得太啰嗦了，你可以使用D=decimal.Decimal; D("1.2") ** D("0.5")。 - recursive
4

首先，您的输入可以大大简化。您不需要读取和解析文件。您只需使用Python表示法声明对象。评估该文件。

b = [
    [4.0, -2.0,  1.0],
    [1.0, +5.0, -3.0],
    [2.0, +2.0, +5.0],
]
y = [ 11.0, -6.0, 7.0 ]

其次，y=-1.2-0.20000000000000001x+0.59999999999999998z并不罕见。对于0.2或0.6，在二进制表示中没有精确的表示方式。因此，显示的值是原始值的十进制近似值而非精确表示。这对于几乎所有类型的浮点处理器都是正确的。

您可以尝试使用Python 2.6的fractions模块。有一个较旧的rational包可能会有所帮助。

是的，将浮点数提高到幂会增加误差。因此，您必须确保避免使用浮点数的最右侧位置，因为那些位大多是噪声。

在显示浮点数时，您必须适当地舍入它们以避免看到噪声位。

>>> a
0.20000000000000001
>>> "%.4f" % (a,)
'0.2000'
我需要从文件中读取输入，因为我的程序必须足够通用，能够处理NxN矩阵。此外，%.4f不是字符串格式化格式吗？它似乎不太干净。 - darudude
1
Python标记法可以描述任意大小的矩阵。为什么要编写自己的解析器？为什么要查找“，”并转换浮点数？只需使用Python符号表示您的矩阵即可。 - S.Lott
所有数字都被转换为字符串以进行打印。您的源代码使用十进制表示法（“0.2”）。Python使用二进制近似值。当您请求输出 - 任何输出 - 它都会被转换回字符串。始终如此。请在该字符串中请求有意义的数字，而不是所有数字。 - S.Lott
4
我建议不要使用decimal模块来处理这样的任务。它的目的更多是处理现实世界中的十进制数字（例如匹配人类记账惯例），有限精度，而不是执行精确的数学运算。在十进制中也有无法准确表示的数字，就像在二进制中一样，并且在十进制中执行算术运算也比替代方案慢得多。
相反，如果你想要精确的结果，应该使用有理数算术。它们将数字表示为分子/分母对，因此可以精确表示所有有理数。如果你只使用乘法和除法（而不是会导致无理数的平方根等操作），你永远不会失去精度。
正如其他人所提到的，Python 2.6将拥有内置的有理数类型，尽管请注意，这不是一个高性能的实现 - 对于速度，最好使用像gmpy这样的库。只需将对float()的调用替换为gmpy.mpq()，你的代码现在应该给出精确的结果（尽管你可能需要将结果格式化为浮点数以供显示）。
下面是稍微整理过的加载矩阵的代码，它将使用gmpy有理数：
def read_matrix(f):
    b,y = [], []
    for line in f:
        bits = line.split(",")
        b.append( map(gmpy.mpq, bits[:-1]) )
        y.append(gmpy.mpq(bits[-1]))
    return b,y
2

我不是回答你的问题，但是与之相关：

#!/usr/bin/env python
from numpy import abs, dot, loadtxt, max
from numpy.linalg import solve

data = loadtxt('gauss.dat', delimiter=',')
a, b = data[:,:-1], data[:,-1:]
x = solve(a, b) # here you may use any method you like instead of `solve`
print(x)
print(max(abs((dot(a, x) - b) / b))) # check solution

例子：

$ cat gauss.dat
4.0, 2.0, 1.0, 11.0
1.0, 5.0, 3.0, 6.0 
2.0, 2.0, 5.0, 7.0

$ python loadtxt_example.py
[[ 2.4]
 [ 0.6]
 [ 0.2]]
0.0
0
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接