浮点数公式的累计误差边界计算

浮点数公式的累计误差边界计算 - 林忠威的博客 | ExtreLin Blog

浮点数公式的累计误差边界计算

<在计算机中很多判断是无法得到一个精确的结果的，这是因为计算底层二进制的形式去描述一个浮点数总是带有误差。在很多应用场景中就需要得到唯一精确的结果，举个简单的例子

double a = 0.2;
double b = 0.3-0.1;
if(a==b)
    std::cout << "a == b" <<std::endl;
else
    std::cout << "a != b" <<std::endl;

运行结果应该是 “a != b”.

在实际应用中线l1、l2和l3相交于一个点p，但是通过l1和l2求得的点p0不一定等于l2和l3求出的点p1，这样就导致后续算法的出现问题。

当然我们可以用一些无限长的精确计算库，如gmp，但是进入精确计算后算法效率就会变得很低。

论文《Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates 》提供了自适应精度的计算的算法。但是论文实在太难了，我只看懂了一部分，所以我只介绍里面最简单的浮点数误差边界的计算。

a ⊕ b = a + b ± err(a ⊕ b )

a ⊙ b = a x b ± err(a ⊙ b )

首先我们计算一个1的的误差边界，其实很简单，当在计算机中 1 + a == 1 则a就是1的误差边界,代码如下：

double get_epsilone()
{
    double epsilon = 1.0;
    double half = 0.5;
    double check = 1.0;
    double lastCheck;
    do{
        lastCheck = check;
        epsilon *=half;
        check = 1.0 + epsilon;
    } while((check != 1.0) &&(lastCheck != check));

    return epsilon;
}

那么上计算a+b的误差边界就是 (abs(a) + abs(b)) * epsilon; 怎么利用这个误差边界呢？

double a,b;

double c = a + b;
double errBound  = (abs(a) + abs(b)) * get_epsilone();
//如果c满足条件则表示c是在误差范围内的
if((c > errBound||(-c > errBound))
    return c;
else
    gmp(a + b); //伪代码，如果不满足就使用gmp去精确计算。

误差边界

结合代码和图，就能明白误差边界的定义，所以我们认为我们计算的数值在两侧误差边界中就表示我们计算的结果可靠。

上面式子是不是太简单了？

我们尝试求解一个简单的几何问题，判断点和直线的位置关系。
判断点和直线的位置关系实际就是求点和直线上的两个点组成的三角形的面积。面积为0表示点在直线上，大于0表示在直线上方，小于0表示在直线的下方。那么判断的式子如下式所示：
$gif.latex?\\det=\left(a_{x}-c_{x}\right)\left(b_{y}-c_{y}\right)-\left(a_{y}-c_{y}\right)\left(b_{x}-c_{x}\right)$