c++ - 哦，我的 OpenMesh vector 算法精度哪里去了？

Question

使用 double s 我希望有大约 15 个小数点的精度。我知道很多十进制数是 not exactly representable in floating point notation , 所以我会得到 1/3 的近似值例如。但是，使用 double我希望得到一个精确到小数点后 15 位的近似值。我还希望在进行算术运算时保持这种准确度。

但是，在下面的示例中，我尝试使用 Heron's formula 来计算三角形的面积和 OpenMesh::Vec3d由 OpenMesh::VectorDataT<double,3> 支持最终得到的结果只能精确到小数点后 5 位。

正确的结果是area = 8.19922e-8 , 但我得到 area=8.1992238711962083e-8 .知道这是从哪里来的吗？

这可能是由 Heron 公式的不稳定性引起的建议是一个很好的建议，但不幸的是，在这个例子中并非如此。我已经为那些可能感兴趣的人添加了计算 Heron 稳定变化的代码。在这个例子中，u.norm()>v.norm()>w.norm() .

#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>

int main()
{
    //triangle vertices
    OpenMesh::Vec3d x(0.051051, 0.057411, 0.001355);
    OpenMesh::Vec3d y(0.050981, 0.057337, -0.000678);
    OpenMesh::Vec3d z(0.050949, 0.057303, 0.0);

    //edge vectors
    OpenMesh::Vec3d u = x-y;
    OpenMesh::Vec3d v = x-z;
    OpenMesh::Vec3d w = y-z;

    //Heron's Formula
    double semiP = (u.norm() + v.norm() + w.norm())/2.0; 
    double area = sqrt(semiP * (semiP - u.norm()) * (semiP - v.norm()) * (semiP - w.norm()) ); 

    //Heron's Formula for small angles
    double areaSmall = sqrt((u.norm() + (v.norm()+w.norm()))*(w.norm()-(u.norm()-v.norm()))*(w.norm()+(u.norm()-v.norm()))*(u.norm()+(v.norm()-w.norm())))/4.0;
}

Answer 1

海伦公式在数值上不稳定。如果你有一个非常“平”的小角度三角形，两条小边的总和几乎就是长边，所以其中一项变得非常小。例如，如果 a 和 b 是小边，

(s - c)

会很小，因为

s = (a + b + c)/2

几乎等于c。

wikipedia article about herons formula提到了一个稳定的替代方案:

排列边使得 a > b > c 并使用

A = 1/4*sqrt((a + (b + c))*(c - (a - b))*(c + (a - b))*(a + (b - c)))