感觉 double 精度不够用啊

V2EX = way to explore

V2EX 是一个关于分享和探索的地方

现在注册

已注册用户请登录

这是一个创建于 431 天前的主题，其中的信息可能已经有所发展或是发生改变。

这几天程序有一个 bug ，查出来是计算一个三角形面积，理论上应该是负数，但是函数算出来是正数，百思不得其解。

后来把计算面积公式里的 double 换成了 doubledouble 类，把精度加倍，结果就正确了。

寻思坐标也不是什么很大的值，怎么 double 就当场挂掉了呢？

=========
测试代码如下:

double buf[6];
buf[0*2+0] = 830776.16442871094; buf[0*2+1] = 430600.56744384766;
buf[1*2+0] = 830776.37707519531; buf[1*2+1] = 430600.68286132813;
buf[2*2+0] = 830776.16430664063; buf[2*2+1] = 430600.56732177734;

int numvertex = 3;

double area = 0.0;
for (int p=numvertex-1,q=0; q<numvertex; p=q++)
{
double a1 = buf[p*2+0] * buf[q*2+1];
double a2 = buf[q*2+0] * buf[p*2+1];
area += a1 - a2;
}

最终 area 结果
如果是 double ，那么值是 0.0000610 （错误）
如果是 doubledouble ，那么值是-0.0000118 （正确）

看起来是计算误差问题，但结果一个正，另一个负，足以把程序逻辑搞得天翻地覆。

23 条回复 • 2023-09-29 20:43:29 +08:00

liuidetmks

2023-09-28 14:16:46 +08:00

可能你需要把工时修改下，尽量减少值相近的数相减，做很小的数的除法

randomxb

2023-09-28 14:17:38 +08:00

浮点数能表示的数在大值上是比较稀疏的, 分布在 0 附近的更集中.
你这两个 double 相乘结果更大了, 所以误差会更大也好解释.
你可以试试把大值规整变小再算下, 比如坐标平移到 0 附近, 结果会更准确.

tool2d

2023-09-28 14:23:13 +08:00

@randomxb 我试了一下把坐标偏移到(0,0)，对这个例子似乎可行。

但总觉得治标不治本，这种精度不够就和地雷差不多，随时会炸。

MicrosoftYK

2023-09-28 14:24:16 +08:00

用高斯面积公式或海伦公式试试，同时可以把坐标值进行归一化，映射到一个较小的范围内，这样计算的时候可能会降低一些舍入误差

lakehylia

2023-09-28 14:25:04 +08:00

如果你需要超高精度的浮点计算，推荐使用专门的超高精度的库。比如 java 里面的 BigDecimal 类。

lakehylia

2023-09-28 14:27:16 +08:00

超高精度的库，内部实现，其实就是以字符来存储每一位十进制数位。然后自己实现一套加减乘除等算法。因为计算机的二进制是没办法精确表达十进制浮点数的。

tool2d

2023-09-28 14:30:41 +08:00

@lakehylia 用高精度类是可以，我就是直觉上认为坐标不算极大数，也不算极小数，用 double 应该能搞定，没想到精度会渗出。

被 double 偷袭了一次，大意了。

MicrosoftYK

2023-09-28 14:31:37 +08:00

同楼上，c/c++可以使用 GMP Library

cy18

2023-09-28 14:42:37 +08:00

这东西不能看绝对精度，要看相对精度。
你这问题相对精度的量级是 0.0000610/(830776.16442871094*430600.56744384766)，太高了。

cy18

2023-09-28 14:45:21 +08:00

0.0000610/(830776.16442871094*430600.56744384766)=1.705 182 951 × 10^(−16),对应倒数在 2^52 跟 2^53 之间，double 的有效位数刚好 2^52 位。

cy18

2023-09-28 14:47:45 +08:00

错了...应当是-0.0000118/(830776.16442871094*430600.56744384766)=3.298 552 355 × 10^(−17)，对应精度 54 位半。

shenjinpeng

2023-09-28 14:49:51 +08:00

建议除了整形(加减乘余) 计算外的所有涉及到小数计算的都用高精度库来处理

MoYi123

2023-09-28 14:59:35 +08:00

感觉代码的逻辑挺奇怪的, 海伦公式求三角面积需要三个点按顺时针(还是逆时针来着?)
否则就是你的例子里的负数, 也就是说调用这个函数前,应该保证三点的顺序.

如果要判断三点的顺序应该用向量求叉积的方式更加自然吧?

nothingistrue

2023-09-28 15:24:35 +08:00

所谓的误差，是计算/测量结果跟预期精确值之间的差异范围。楼主的计算有点高深，没看出来是干啥的。但大概他是要计算 a1 比 a2 下。那么当你考虑误差以后，结果应该是 a1 - a2 的绝对值，小于误差范围，而不该是 a1 - a2 < 0 ，后者是精确值，不是误差范围。

浮点数的精度仅代表当前值的误差范围，而误差会随着不同的计算逻辑而发生变化，所以浮点计算最后的总误差，是要根据不同的算法做不同的分析的。没有人会喜欢这么干，所以如果需要精确计算，不要用浮点数，要用整数或者特定的十进制数据。

Masoud2023

2023-09-28 15:40:49 +08:00

什么语言？

我的建议是涉及到这种精密计算有条件都算专门的计算库（ BigDecimal ？）

能减少很多坑。

msg7086

2023-09-28 16:07:02 +08:00

这是相当于在两个奥运场馆大小的区域用游标卡尺比大小啊，这么精细的操作 double 确实是不够的。

misdake

2023-09-28 16:09:03 +08:00

不去改计算方法来提高精度的话，接下来 doubledouble 也会不够用的。

williamscorn

2023-09-28 16:10:02 +08:00

参考
https://en.wikipedia.org/wiki/Double-precision_floating-point_format
可知，double 一般只有 15~17 位的有效十进制位，测试用例计算两点叉积中间过程如下：

double
357732783707.917969-357732779387.523071 = 4320.394897
357732779286.109924-357732783655.354370 = -4369.244446
357732687769.262695-357732687720.413086 = 48.849609

long double
357732783707.917946-357732779387.523051 = 4320.394895
357732779286.109919-357732783655.354386 = -4369.244467
357732687769.262668-357732687720.413108 = 48.849560

可以看到在 16 位之后，double 和 long double 的结果已经不同了，这导致了最后结果的不同。

PS：这样计算完，还需要再*0.5 再取绝对值，才是三角形面积吧。

codeTempo

2023-09-28 16:21:59 +08:00 via Android

两个相近的浮点数做差会有精度上的问题
https://en.m.wikipedia.org/wiki/Catastrophic_cancellation

geelaw

2023-09-28 16:40:56 +08:00 via iPhone

第一个考虑的方向是条件数，从三个平直系坐标计算三角形面积，在所给的输入处的条件数是 400 亿，这说明计算结果只能期待 <6 位有效数字。

第二个考虑的方向是误差累积，楼主的算法，绝对误差是 (sum |x_i| + sum |y_i|) * eps （这个界比较松）。

如果先把任意一个点平移到原点，则绝对误差是 (sum |x_i - x_1| + sum |y_i - y_1|) * 2eps ，是更好的策略。

zmxnv123

2023-09-28 16:44:26 +08:00

为什么三角形的面积可以是负数

MeteorCat

2023-09-29 10:15:14 +08:00 via Android

用 bigdecimal 相关库绝对没错，其他技巧都是假的

realpg

2023-09-29 20:43:29 +08:00

不是科班出身就容易踩这种坑
能不用 float 就不用 float
必须得用 float ，每次任意运算之前都得先问问自己精度的问题