除双精度数时出现意外的精度损失

6

我有一个函数getSlope,它需要4个double类型的参数,并返回另一个double类型的值,该值是根据给定的参数以以下方式计算得出的:

double QSweep::getSlope(double a, double b, double c, double d){
double slope;
slope=(d-b)/(c-a);
return slope;
}

问题在于调用此函数时,例如带有参数:
getSlope(2.71156, -1.64161, 2.70413, -1.72219);

返回的结果是:

10.8557

这对我的计算结果来说并不理想。 我使用Mathematica计算出相同参数的斜率结果为:

10.8452

或者使用更多数字来提高精度:

10.845222072678331.

我的程序返回的结果在进一步计算中不够好。此外,我不明白程序如何从10.845222072678331(假设这是除法的近似结果)开始返回10.8557。我该如何获得正确的除法结果?
提前感谢您, Madalina
我使用命令行打印结果:
std::cout<<slope<<endl;

也许我的参数可能不够好,因为我是从另一个程序中读取它们的(该程序计算图形;在我读取这些参数后,我只是将它们显示出来以查看其值,但可能显示的向量与计算出的值的内部精度不同。我不知道,这真的很奇怪。会出现一些数值误差..)
当计算我要读取参数的图形时,使用了一些用C++编写的数值库(带有模板)。此计算过程中未使用OpenGL。
谢谢, madalina

检查方法在汇编中是如何编译的。我相信你可以在调试器中完成这个操作(至少在Visual Studio中)。 - Saulius Žemaitaitis
我的Windows计算器给出的结果和Mathematica一样好 :D - klew
哪个编译器,哪个平台? - peterchen
8个回答

7
我已经尝试了使用float而不是double,结果为10.845110。与Madalina的结果相比仍然更好。
编辑:
我想我知道为什么会有这样的结果。如果您从其他地方获取a、b、c和d参数并打印它们,它会给您舍入后的值。然后,如果您将其放入Mathematica(或calc ;))中,它将给出不同的结果。
我尝试稍微更改了其中一个参数。当我这样做时:
double c = 2.7041304;

我得到了10.845806。我只是将0.0000004加入c!因此,我认为你的“错误”并不是错误。使用更好的精度打印a、b、c和d,然后将它们放入Mathematica中。

5

你的项目中是否使用了DirectX或OpenGL?如果是,它们可能会关闭双精度,导致出现奇怪的结果。

你可以通过以下方式检查你的精度设置:

std::sqrt(x) * std::sqrt(x)

结果必须非常接近于x。 我很久以前遇到了这个问题,花了一个月的时间检查所有公式。但是后来我找到了。
D3DCREATE_FPU_PRESERVE

他们究竟是如何做到的呢? - anon
在初始化Direct3D时有一些选项。我不记得名字了,但我花了一个月的时间检查所有文凭公式,只有当我用"sqrt(x)*sqrt(x)"进行简单检查时,精度才会大大降低,除非我关闭该选项。 - Mykola Golubyev
当使用VS2008中的标准Win32控制台应用程序编译时,它会给出正确的答案。我同意并认为这是编译器设置的问题。 - Binary Worrier

5
以下代码:
#include <iostream>
using namespace std;

double getSlope(double a, double b, double c, double d){
    double slope;
    slope=(d-b)/(c-a);
    return slope;
}

int main( ) {
    double s = getSlope(2.71156, -1.64161, 2.70413, -1.72219);
    cout << s << endl;
}

使用g++编译器,结果为10.8452。请问您在代码中如何打印输出结果?


无论你如何打印10.845222072678331,它都不会四舍五入或截断为10.8557。 - Pete Kirkham

3
这里的问题是(c-a)很小,因此在浮点运算中固有的舍入误差在这个例子中被放大了。一般的解决方案是重新设计你的方程,以便你不需要除以一个小数,但我不确定你在这里该如何做。
编辑:Neil在他对这个问题的评论中是正确的,我使用双精度浮点数在VB中计算得到与Mathematica相同的答案。

2
您得到的结果与32位算术一致。如果不知道更多关于您的环境信息,就无法建议应该做什么。
假设所示代码是正在运行的代码,即您没有将任何内容转换为字符串或浮点数,则在C++中没有修复方法。它在您所展示的代码外部,并取决于环境。
由于Patrick McDonald和Treb都提出了输入精度和a-c误差的准确性,我想看一下这个问题。查看舍入误差的一种技术是区间算术,它使值表示的上限和下限显式(它们在浮点数中是隐含的,并且固定为表示的精度)。通过将每个值视为上限和下限,并通过扩展表示中的误差的边界(对于双倍精度值x,大约为x * 2 ^ -53),您会得到一个结果,该结果给出了一个值的下限和上限,考虑到最坏情况下的精度误差。
例如,如果您有一个范围在[1.0,2.0]中的值,并从其中减去一个范围在[0.0,1.0]中的值,则结果必须位于[below(0.0),above(2.0)]范围内,因为最小结果为1.0-1.0,而最大结果为2.0-0.0。below和above相当于floor和ceiling,但是对于下一个可表示的值而不是整数。
使用表示最坏情况下双重舍入的区间:
getSlope(
 a = [2.7115599999999995262:2.7115600000000004144], 
 b = [-1.6416099999999997916:-1.6416100000000002357], 
 c = [2.7041299999999997006:2.7041300000000005888], 
 d = [-1.7221899999999998876:-1.7221900000000003317])
(d-b) = [-0.080580000000000526206:-0.080579999999999665783]
(c-a) = [-0.0074300000000007129439:-0.0074299999999989383218]

to double precision [10.845222072677243474:10.845222072679954195]

因此,尽管 c-a 与 c 或 a 相比较小,但仍然大于双重舍入,因此如果您使用的是最糟糕的双精度舍入,则可以相信该值精确到12个数字-10.8452220727。 您已经失去了一些双精度数字,但仍然可以更加准确地处理您的输入数据。
但是,如果输入只有几个有效数字,而不是双精度值2.71156 +/- eps,则输入范围将为[2.711555, 2.711565],因此您将得到以下结果:
getSlope(
 a = [2.711555:2.711565], 
 b = [-1.641615:-1.641605], 
 c = [2.704125:2.704135], 
 d = [-1.722195:-1.722185])
(d-b) = [-0.08059:-0.08057]
(c-a) = [-0.00744:-0.00742]

to specified accuracy [10.82930108:10.86118598]

这是一个更广泛的范围。

但您需要费心去跟踪计算的准确性,而浮点数中固有的舍入误差在此示例中并不重要 - 即使采用最坏情况下的双精度舍入,它也精确到12个数字。

另一方面,如果您的输入只知道6个数字,那么无论您得到10.8557还是10.8452都没有关系。两者都在[10.82930108:10.86118598]之内。


1

最好也打印出参数。我猜你是用十进制表示法传递参数,但每个参数都会失去精度。问题在于1/5在二进制中是一个无限序列,所以例如0.2变成了.001001001.... 此外,在将二进制浮点数转换为十进制文本表示时,小数部分会被截断。

此外,有时编译器会选择速度而不是精度。这应该是一个记录在案的编译器开关。


0

Patrick 似乎是对的,认为 (c-a) 是主要原因:

d-b = -1.72219 - (-1.64161) = -0.08058

c-a = 2.70413 - 2.71156 = -0.00743

S = (d-b)/(c-a)= -0.08058 / -0.00743 = 10.845222

你开始使用六位数精度,通过减法运算使其降至 3 和四位数。我最好的猜测是,因为数字 -0.00743 不能在双精度中精确表示,所以你失去了额外的精度。尝试使用具有更大精度的中间变量,像这样:

double QSweep::getSlope(double a, double b, double c, double d)
{
    double slope;
    long double temp1, temp2;

    temp1 = (d-b);
    temp2 = (c-a);
    slope = temp1/temp2;

    return slope;
}

你似乎混淆了精度(数字的表示方式)和准确性(值的公差)。在C++中,无论你将double指定为2.70413还是2.7041300000,对结果没有影响。 - Pete Kirkham
@Pete Kirkham:在double中无法精确表示例如0.1的值,因此将其存储在范围更大的变量中可能会产生不同的结果。 - Treb
正如我在我的回答中所展示的,“您从六位数字精度开始,”与OP代码给出的结果无关。它与您应该关心多少个结果数字有关,但是如果使用64位双精度计算结果,则(c-a)不是错误的原因。 - Pete Kirkham

-1

虽然学术讨论对于了解编程语言的限制非常有帮助,但您可能会发现解决问题最简单的方法是使用任意精度算术的数据结构。

这将会增加一些额外开销,但您应该能够找到某些准确性相当可靠的解决方案。


1
推荐使用任意精度算术,虽然在StackOverflow上很受欢迎,但并不是每个关于浮点数计算的问题的最佳答案。 - quant_dev
这是正确的,但通常最简单可行的解决方案仍然是最好的。虽然有时候有更好、更快和更复杂的方法来完成任务,但即使只有简单性在软件项目中也有很大的价值。 - Ian Gilham

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接