double - 及其用法

评分：3.1/5 (8 票)

本文介绍 "double" 算术的一些特性，以及你可以和不可以做什么。
注意：需要 UTF-8 支持数学符号。

0. 动机
=============
在进行浮点运算时，可能会出现一些奇怪的效果，例如

for(double i=0.0; i!=5.5; i+=0.5)
{
        ... 
}

由于舍入误差，可能永远不会终止。特别是测试 (不) 等值通常不会得到期望的结果。一个常见的解决方法是使用一个 epsilon

bool is_equal(double d1, double d2) 
{
        if(abs(d1-d2)<epsilon)
                return true;
        return false;
}

然而，这个 "解决方案" 只在少数情况下适用 - 这样定义的等价性不再是等价关系（失去了传递性：a-b<ε 且 b-c<ε 并不能推导出 a-c<ε - 哎呀）。另一个问题是：epsilon 应该有多大？

计算机应该计算出正确的结果，或者说明它的不确定性。使用 epsilon 通常无法实现这一点。这就是为什么你应该对本文的其余部分感兴趣 ;-)

1. 数值类型
==================
IEEE 754 double 值由以下几部分组成：
1 位符号 (+ 或 -)
11 位指数（范围：0 到 1047，0 和 1047 是 "保留" 的）
52 位有效数字（也称为 "尾数"）

则表示的数字为

[符号] 2^{指数-偏差[=1023]}*1[有效数字]

其中 1 是一个隐含位，在有效数字位开始之前设置为 1。

对于归一化数字，有效数字的第一位是 1 且不保存，因此有效数字扩展到 53 位（显然，前导零可以省略，因此数字以非零开头，在二进制系统中碰巧总是 1）。

当且仅当 0 < 指数 < 2047 时，该值才被归一化。如果指数为零，IEEE 使用 "非归一化" 数字，则值为 [符号]2^-1022*[有效数字]，但不包含有效数字中的前导 1。请注意，在非归一化数字范围内的计算通常非常慢（指数中的 2047 保留给 "NaN (非数字)、±inf 及类似项"）。

如果你出于某种原因需要位访问，可以在 *大端序* 机器上使用

union ieee_double
{
   double   a;
   struct { unsigned s : 1;
            unsigned e :11;
            unsigned h :20;
            unsigned l :32;
          } b;
};

小端序需要单独处理，但方式类似。

2. 操作
=============
假设舍入设置为 "舍入到最近" 和 "偶数舍入断点"（默认）

设 ∘ 表示 {+,-,*,/} 中的一个运算，设 ⊚ 表示相应的机器运算，即 {⊕,⊖,⊛,⊘} 中的一个。

设 x、y 和 z 为 IEEE 754 double 值。
设 a 为实数。
设 fl(a) 为最接近 a 的 IEEE 754 double 值。

IEEE 754 现在保证
(a) fl(x ∘ y) = x ⊚ y
推论
x ⊖ y = 0 ⇔ x = y
x ⊕ y = y ⊕ x
x ⊛ y = y ⊛ x
(b) fl(sqrt(x)) = msqrt(x)，其中 msqrt 是 "机器 sqrt" 操作。
此外，还定义了一些 "非数字" 结果，例如 NaN（非数字）和 +/- 无穷大。

然而，在一般情况下，由于舍入误差
a ⊚ ( b ⊚ c ) ≠ ( a ⊚ b ) ⊚ c

3. Sterbenz 引理
=================
在进行 double 运算时，你最重要的朋友之一是 Sterbenz 引理，它可以表述为：
如果 y/2 ≤ x ≤ y 则
x ⊖ y = x - y
如果 x ⊖ y 不下溢

或者，用 "普通英语" 表述：如果 y 在 x/2 和 x 之间，那么机器减法 x ⊖ y 就会给你精确的结果。

4. 误差界
===============
给定 (1) 和 (2) 中的性质以及由此产生的引理 (3)，可以为使用 +,-,*,/ 和 'sqrt' 的任何公式计算误差界。然而，这种计算本身很繁琐且容易出错，因此通常使用近似值。

利用 Sterbenz 引理，可以证明一些更简单的多项式求值误差界。我不会在这里展示它们，但你可以搜索例如：
- Fortune-van-Wyk 界
- Mehlhorn-Näher 界
- Mehlhorn-Neumaier 界
来找到它们。（注意：据我所知，所有这些论文都仅证明了它们对于整数操作数（在 doubles 中表示）的正确性。然而，它们对任何 double 输入都正确。）所有这些公式一个值得注意的美妙之处在于，实际的误差界可以使用 double 算术来计算。

5. 浮点过滤器
=========================
记住动机示例：等值测试。给定两个多项式 P 和 Q，现在可以在某些情况下确定 P(x) = Q(x)。

// equality test:
if(abs(P(x1) - Q(x2)) > ErrorBound(max(x1,x2)))
        return false; /* we know that for sure */
else
        // ???

（请注意，如果 P 或 Q 是多元的，max 最好使用 x1 和 x2 的上确界范数...）

如果你在编译时（或程序启动时）知道操作数的大小并且可以预先计算它，那么在这种情况下，你只需要进行一次简单的测试即可获得精确性。
问题：'???' 是什么？有几种选择。你可以使用任意精度数字类型的精确算术，你可以使用浮点展开，或者你可以执行惰性自适应求值（请注意，速度从左到右增加，但实现难度（以及想出正确公式的难度）也随之增加）。
这些方法超出了本文的范围（如果时间允许，我可能会写一篇续篇 ;-））。
在某些领域，另一种可能性是使用受控扰动，这意味着输入被转换成不会出现等值且可以显示不等值的情况（不精确，但可能鲁棒且足够，因为我们保证了某种属性（即，我们不仅仅说 "不相等"，我们还确保它不是）。

6. 结论
=============
虽然计算误差界需要额外的工作，但运行时效率是一个小问题（在大多数情况下，只需一次比较即可证明结果的正确性）。同时，由于错误地假设等值而导致的许多退化情况可以避免，通常可以显著提高性能。

________________________________________________________________________________
如果您有任何评论、建议或改进，或者发现了拼写错误（或多个），请告知我。

C++

文章

double - 及其用法