这是一个创建于 212 天前的主题，其中的信息可能已经有所发展或是发生改变。
这几天程序有一个 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 （正确）
看起来是计算误差问题，但结果一个正，另一个负，足以把程序逻辑搞得天翻地覆。
后来把计算面积公式里的 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 （正确）
看起来是计算误差问题，但结果一个正，另一个负，足以把程序逻辑搞得天翻地覆。
 |
1
liuidetmks 212 天前
可能你需要把工时修改下，尽量减少 值相近的数相减，做很小的数的除法
|
 |
2
randomxb 212 天前  2
浮点数能表示的数在大值上是比较稀疏的, 分布在 0 附近的更集中.
你这两个 double 相乘结果更大了, 所以误差会更大也好解释. 你可以试试把大值规整变小再算下, 比如坐标平移到 0 附近, 结果会更准确. |
 |
4
MicrosoftYK 212 天前
用高斯面积公式或海伦公式试试，同时可以把坐标值进行归一化，映射到一个较小的范围内，这样计算的时候可能会降低一些舍入误差
|
 |
5
lakehylia 212 天前 1
如果你需要超高精度的浮点计算，推荐使用专门的超高精度的库。比如 java 里面的 BigDecimal 类。
|
|
6
lakehylia 212 天前
超高精度的库，内部实现，其实就是以字符来存储每一位十进制数位。然后自己实现一套加减乘除等算法。因为计算机的二进制是没办法精确表达十进制浮点数的。
|
 |
7
tool2d OP |
|
8
MicrosoftYK 212 天前
同楼上，c/c++可以使用 GMP Library
|
 |
9
cy18 212 天前
这东西不能看绝对精度，要看相对精度。
你这问题相对精度的量级是 0.0000610/(830776.16442871094*430600.56744384766)，太高了。 |
|
10
cy18 212 天前
0.0000610/(830776.16442871094*430600.56744384766)=1.705 182 951 × 10^(−16),对应倒数在 2^52 跟 2^53 之间，double 的有效位数刚好 2^52 位。
|
|
11
cy18 212 天前 1
错了...应当是-0.0000118/(830776.16442871094*430600.56744384766)=3.298 552 355 × 10^(−17)，对应精度 54 位半。
|
 |
12
shenjinpeng 212 天前
建议除了 整形(加减乘余) 计算外的所有涉及到小数计算的都用高精度库来处理
|
 |
13
MoYi123 212 天前
感觉代码的逻辑挺奇怪的, 海伦公式求三角面积需要三个点按顺时针(还是逆时针来着?)
否则就是你的例子里的负数, 也就是说调用这个函数前,应该保证三点的顺序. 如果要判断三点的顺序应该用向量求叉积的方式更加自然吧? |
 |
14
nothingistrue 212 天前
所谓的误差，是计算/测量结果跟预期精确值之间的差异范围。楼主的计算有点高深，没看出来是干啥的。但大概他是要计算 a1 比 a2 下。那么当你考虑误差以后，结果应该是 a1 - a2 的绝对值，小于误差范围，而不该是 a1 - a2 < 0 ，后者是精确值，不是误差范围。
浮点数的精度仅代表当前值的误差范围，而误差会随着不同的计算逻辑而发生变化，所以浮点计算最后的总误差，是要根据不同的算法做不同的分析的。没有人会喜欢这么干，所以如果需要精确计算，不要用浮点数，要用整数或者特定的十进制数据。 |
 |
15
Masoud2023 212 天前
什么语言？
我的建议是涉及到这种精密计算有条件都算专门的计算库（ BigDecimal ？） 能减少很多坑。 |
 |
16
msg7086 212 天前
这是相当于在两个奥运场馆大小的区域用游标卡尺比大小啊，这么精细的操作 double 确实是不够的。
|
 |
17
misdake 212 天前
不去改计算方法来提高精度的话，接下来 doubledouble 也会不够用的。
|
 |
18
williamscorn 212 天前 1
参考
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 再取绝对值，才是三角形面积吧。 |
 |
19
codeTempo 212 天前 via Android
两个相近的浮点数做差会有精度上的问题
https://en.m.wikipedia.org/wiki/Catastrophic_cancellation |
 |
20
geelaw 212 天前 via iPhone
第一个考虑的方向是条件数，从三个平直系坐标计算三角形面积，在所给的输入处的条件数是 400 亿，这说明计算结果只能期待 <6 位有效数字。
第二个考虑的方向是误差累积，楼主的算法，绝对误差是 (sum |x_i| + sum |y_i|) * eps （这个界比较松）。 如果先把任意一个点平移到原点，则绝对误差是 (sum |x_i - x_1| + sum |y_i - y_1|) * 2eps ，是更好的策略。 |
 |
21
zmxnv123 212 天前
为什么三角形的面积可以是负数
|
 |
22
MeteorCat 211 天前 via Android
用 bigdecimal 相关库绝对没错，其他技巧都是假的
|
 |
23
realpg 211 天前
不是科班出身就容易踩这种坑
能不用 float 就不用 float 必须得用 float ，每次任意运算之前都得先问问自己精度的问题 |
Select Language