舍入误差

硬件浮点数的四舍五入方式异常简单：只有在操作的结果不能精确表示时，四舍五入才会发生，而且默认情况下，会采取四舍五入到最接近的可表示数字（在平局的情况下，优先选择以零结尾的数字）。 考虑以下代码片段：

float x = 0;
for (int i = 0; i < (1 << 25); i++)
    x++;
printf("%f\n", x);

结果输出1677216=2241677216 = 2^{24}，而不是数学上期望的结果225=335544322^{25} = 33554432。为什么呢？

当重复递增一个浮点数xx时，我们最终会到达这样一个状态，浮点数变得如此之大以至于(x+1)(x+1)会被舍入回xx。递增得到的第一个这样的数是2242^{24}（即尾数位数加1），因为2242^{24}作为浮点数递增（尾数末位加1）的下一个值是：

(float)224+1=224⋅1.0...01⏟×23=1677218(float)2^{24} + 1 = 2^{24} \cdot 1.\underbrace{0...01}_{×23} = 1677218

即用32位浮点数无法精确表示1677217这个整数，而由于1677217到2242^{24}和224+12^{24}+1的距离相同，所以根据之前提到的“平局断定规则”，1677217会被向下舍入到2242^{24}。同时，任何小于2^24的数字，它们的增量可以被精确地表示出来，也就是说，在这个范围内的数在计算过程中不需要进行四舍五入的操作。

舍入误差和操作顺序

浮点计算的结果可能会因为操作顺序的不同而不同，尽管在纯数学意义上加法和乘法运算是可交换和可结合的，但它们的舍入误差却不是。

例如，假设我们有三个浮点变量xx、yy和zz，(x+y+z)(x+y+z)的结果会因为加法的顺序不同而不同。这种不可交换性原则适用于其他大多数浮点运算。

编译器不能产生不符合规范的结果，所以这个烦人的细微差别禁止了一些涉及到在算术中重新排列操作数的可能优化。在GCC和Clang中，你可以用-ffast-math标志来禁用严格遵守规范。如果我们添加这个标志并重新编译上面的代码片段，它运行得会快很多，并且也会输出正确的结果33554432（但你需要注意，编译器也可能选择一条计算精度较低的路径）。

舍入模型

除了默认的四舍五入方法（也被称为银行家舍入（Banker’s rounding））外，你还可以设置其他4种舍入逻辑模式：

  • “四舍五入，完美平衡时总是远离0舍入”：这是最常用的舍入模式，当某个数的小数部分刚好处于中间值时（即0.5），按照远离0的方向进行舍入。
  • “向上舍入（朝向+∞+∞，负数结果因此向零进行舍入）”：这种模式会将所有数都向上舍入，如果是正数则向大的方向舍入，如果是负数则更接近0，也就是向大的方向舍入。
  • “向下舍入（朝向−∞-∞，负数结果因此远离零进行舍入）”：这个模式的操作刚好与向上舍入相反，所有数都向下舍入，如果是正数则更接近0，如果是负数则向小的方向舍入。
  • “向零舍入（二进制结果的截断）”：这种模式实际上等于对数值进行截断，去掉小数部分，不考虑是否需要舍入。例如，3.7和-3.7在这种模式下都将被舍入为3和-3，而不管小数部分的大小。

举个例子，如果在运行上述循环之前调用fesetround(FE_UPWARD)这个函数（也就是设置为向上舍入模式），输出的结果不是2242^{24}，也不是2252^{25}，而是67108864，也就是2262^{26}。 这是因为当我们计算到2242^{24}时，(x+1)(x+1)开始向上舍入到最接近的可表示数值(x+2)(x+2)，以至于我们在原本一半的迭代次数内就到达了2的25次方。接着，(x+1)(x+1)向上舍入到(x+4)(x+4)，我们的计算速度变得更快，一次等于计算四个单位。

切换舍入模式的其中一个应用是用来诊断数值稳定性。如果一个算法的结果在切换舍入到正无穷大和负无穷大之间时有大幅度变化，那么这意味着这个算法可能对舍入误差有较高的敏感性。

这种测试方法通常比将所有的计算切换到更低精度并检查结果是否有大变化要好，因为默认的“向最接近的数舍入”策略会在足够的平均后收敛到正确的“期望”值：在一半的迭代次数中，误差会向上舍入，而在另一半迭代中，误差会向下舍入——所以从统计上看，这些误差会相互抵消。

测量误差

尽管执行复杂的运算任务（例如自然对数、平方根）想要保证精度看似困难，但事实是硬件能保证提供最高的计算精度，结果准确度极高。这种特性使得分析舍入误差变得异常简单。

有两种衡量计算错误的方法：

  • 创建硬件或符合规范的精确软件的工程师关心的是最后一位的单位（units in the last place，ulps），即两个数字之间的距离，即精确实数和实际计算结果之间可以容纳多少可表示的数字。
  • 研究数值算法的人则关心相对精度，这是由误差绝对值除以真实结果得出的数值：∣v−v,v∣|\frac{v-v^,}{v}|。

无论哪种情况，一种常见的策略是假设最坏的情况，并给出其上限。

也就是说，在执行单个基础算术操作时，最糟糕的情况就是结果舍入到最近的可表示数。这意味着误差不超过0.5个单位最后一位（ulps）。至于相对误差，我们可以定义一个名为机器精度（machine epsilon）的数字ϵ\epsilon，为1和下一个可被计算机表示的数的差值，这个差值等于2的负幂次数，而这个负幂次数等于为浮点数的尾数（mantissa）所分配的位数。

这意味着，在进行一次算术操作后，你得到了某个结果xx，真正的结果应该在一个范围内：

[x⋅(1−ϵ), x⋅(1+ϵ)][x \cdot (1 - \epsilon),\space x\cdot (1 + \epsilon)]

在进行浮点数计算的结果分析时，必须要知道的一点是，这种计算总是存在一定的误差。当你需要根据这些结果做出“是或否”的决策时，这个因素就显得尤为重要。举个例子，你可能需要通过浮点计算检查两个值是否相等：

const float eps = std::numeric_limits<float>::epsilon; // ~2^(-23)
bool eq(float a, float b) {
    return abs(a - b) <= eps;
}

eps的值应根据应用的需要来确定。不同的应用情境对精度的要求可能会有所不同。 上面的eps值（对应float的机器精度）适用于最多一次的浮点运算。这是因为每次进行浮点运算，尤其是进行多次连续浮点运算时，都可能带来新的误差。因此，如果你进行了多次的浮点运算，使用上述的机器精度可能就不再适用，你可能需要选择一个更大的eps值来反映运算中可能产生的更大误差。

区间算术（Interval Arithmetic）

一个算法如果其在计算过程中的误差不会显著增加，无论误差的来源是什么，那么这个算法就可以被称为数值稳定（numerical stable）的。只有在问题本身条件良好（well-conditioned）的情况下，也就是说，如果输入数据发生微小变化，解决方案只会发生微小变化的情况（即对输入的微小变化不敏感）下，这才可能发生。

在分析数值算法时，采用与实验物理相同的方法往往非常有用：我们不是处理未知的实数值，而是处理它们可能所在的区间。

例如，考虑一个连续操作链，我们在其中连续地将一个变量乘以任意实数：

float x = 1;
for (int i = 0; i < n; i++)
    x *= a[i];

在第一步乘法之后，实际的乘积应该是xx，但是由于计算误差，我们得到的计算结果是x⋅(1+ϵ)x\cdot(1+ϵ)。 接下来，我们每做一个额外的乘法操作，这个误差的上限就要乘以另一个(1+ϵ)(1+ϵ)。通过归纳法，可以看出在n次连续的乘法之后，计算出来的值将在(1+ϵ)n=(1+nϵ+O(ϵ2))(1 + ϵ)^n = (1+nϵ+O(ϵ^2))这个范围内上下浮动。这个等式是基于泰勒级数展开得来的。O(ϵ2)O(ϵ^2)表示的是ϵϵ的平方级别的项。当ϵϵ趋近于0，这个项的值也会趋近于0，所以在ϵϵ非常小的情况下，这个项对结果的影响可以不考虑。

这意味着相对误差为O(nϵ)O(n\epsilon)，这一般是没问题的，因为通常n≪1ϵn \ll \frac{1}{\epsilon}。

举个数值不稳定的例子，考虑如下函数：

f(x,y)=x2−y2f(x,y) = x^2 - y^2

假设x>yx > y，这个函数返回的最大值约为：

x2⋅(1+ϵ)−y2⋅(1−ϵ)x^2\cdot(1 + \epsilon) - y^2\cdot(1-\epsilon)

对应的绝对误差：

x2⋅(1+ϵ)−y2⋅(1−ϵ)−(x2−y2)=(x2+y2)⋅ϵx^2\cdot(1 + \epsilon) - y^2\cdot(1-\epsilon) - (x^2 - y^2) = (x^2 + y^2) \cdot \epsilon

因此相对误差为：

x2+y2x2−y2⋅ϵ\frac{x^2 + y^2}{x^2 - y^2}\cdot \epsilon

如果xx和yy在量级上接近，则误差为O(ϵ⋅∣x∣)O(\epsilon \cdot |x|)。

在直接计算的情况下，减法“放大”了平方的误差。但这可以通过使用以下公式来解决：

f(x,y)=x2−y2=(x+y)⋅(x−y)f(x,y) = x^2 - y^2 = (x + y)\cdot(x - y)

这里很容易得到误差上限为ϵ⋅∣x−y∣\epsilon \cdot |x - y|。同时这样计算也更快，因为这样需要2次加法和一次乘法：比原来多了一次快速的加法操作，少了一次慢速的乘法操作。

Kahan累加算法（Kahan Summation）

进行一长串的计算操作本身并不会引发问题，但是在处理大规模和小规模数值的加减运算时就会出现问题。因为大数和小数在进行加减运算时，往往会导致精度损失。因此，处理这类问题的通用方法是尽量使大的数与大的数进行运算，小的数与小的数进行运算。

考虑标准的累加算法：

float s = 0;
for (int i = 0; i < n; i++)
    s += a[i];

因为我们在进行求和运算，而非乘法运算，因此，它的相对误差不再仅仅受到O(ϵ⋅n)O(ϵ⋅n)（即，误差项ϵϵ和求和项数nn的乘积）的上限约束。反而，在很大程度上，它取决于输入的数据。

考虑最极端的一种情况，如果第一个值为2242^{24}，其他值为1，则不管nn为多少，累加和都为2242^{24}，可以通过执行以下代码进行验证，你会发现两次打印的结果均为16777216=22416777216=2^{24}：

const int n = (1<<24);
printf("%d\n", n);

float s = n;
for (int i = 0; i < n; i++)
    s += 1.0;

printf("%f\n", s);

float浮点数只有23位尾数位（mantissa bits），意味着224+12^{24}+1是其第一个不能精确表示，必须向下取整的整数，因此每当我们试图给s=224s = 2 ^ {24} 加11，就会产生舍入误差。 如果最后一个数恰好为−224-2^{24}，则误差就会接近无限大。

一种显而易见的解决方案是将变量的类型转换为更大的类型，如double，但这个方法的可扩展性不强。一种更优雅的解决方案是将无法添加到原变量的部分另行存储在一种变量中，然后将这个变量添加到下一个变量。

float s = 0, c = 0;
for (int i = 0; i < n; i++) {
    float y = a[i] - c; // c is zero on the first iteration
    float t = s + y;    // s may be big and y may be small, losing low-order bits of y
    c = (t - s) - y;    // (t - s) cancels high-order part of y
    s = t;
}

这个技巧被称为Kahan求和。它的相对误差上限为2ϵ+O(nϵ2)2\epsilon + O(n\epsilon^2)：第一项来自最后一次求和，第二项是因为我们在每次累加上的误差都小于ϵ\epsilon。

当然，常用的处理计算精度问题的方法是使用更高精度的数据类型，比如double型，这也相当于将机器精度提高到平方级别。还可以通过某种方式将两个double型变量捆绑在一起：一个存储数值，另一个存储不可表示的误差，这样两个变量就可以表示数值a+ba+b。这种方法被称为“双双精度”（double-double）运算。此外，这种方法还可以推广到“四倍双精度”（quad-double）和更高精度的运算。