在嵌入式系统中进行更理想的浮点计算的技巧 理工大学, 新加坡, 嵌入式, 表达式, 技巧

pengpengpang

浮点算法不遵循整数算法规则，但利用 FPGA 或者基于 FPGA 的嵌入式处理器不难设计出精确的浮点系统。
作者：Sharad Sinha
在读博士
新加坡南洋理工大学
sharad_sinha@pmail.ntu.edu.sg
工程人员一看到浮点运算就会头疼，因为浮点运算用软件实现速度慢，用硬件实现则占用资源多。理解和领会浮点数的最佳方法是将它们视为实数的近似值。实际上这也是开发浮点表达式的目的。正如老话所说，真实的世界是模拟的，有许多电子系统在与真实世界的信号打交道。因此，对于设计这些系统的工程师来说，有必要理解浮点表达式以及浮点计算的优势和局限性。这将有助于我们设计可靠性更高、容错性更好的系统。
首先深入了解一下浮点运算。通过一些示例计算就可以看到浮点运算不如整数运算直接，工程师在设计使用浮点数据的系统时必须考虑这一点。这里有个重要的技巧：用对数来运算极小的浮点数据。我们的目的是熟悉一些数值运算的特点，把重点放在设计问题上。本文结尾列出的参考文献中有更加深入的介绍。
就设计实现而言， 设计人员可以使用赛灵思LogiCORE™ IP Floating-Point Operator 生成和实例化RTL 设计中的浮点运算符。如果采用 FPGA 中的嵌入式处理器构建浮点软件系统，则可以使用 Virtex-5 FXT 中面向PowerPC 440 嵌入式处理器的 LogicCORE IP Virtex®-5辅助处理器单元 (APU) 浮点单元 (FPU)。
浮点运算 101
IEEE 754-2008 是现行的浮点运算 IEEE 标准。(1)它取代了该标准的 IEEE 754-1985 和 IEEE 754-1987 版本。该标准中的浮点数据表达规则如下：
 要求零带符号：+0，-0
 非零浮点数可表达为 (-1)^s x b^e x m，这里：s 为 +1 或者 -1，用于表明该数为正数还是负数，b是底数（2 或者 10），e 是指数，m 则是以 d0.d1d2—dp-1 形式表达的数，这里 di 对以 2 为底数的情况可以是 0 或者 1，对以 10 为底数的情况可以是0 和 9 之间的任意值。请注意小数点应紧跟 d0。
 分正负的极限：+∞，-∞。
 非数值，有两种形式：qNaN（静态）和 sNaN（信号）。
表达式 d0.d1d2—dp-1 指“有效值”，e 为“指数”。有效值总共有 p 位数，p 即为表达式的精度。IEEE754-2008 定义了五种基本的表达式格式，三种用于 2 为底数的情况，两种用于 10 为底数的情况。标准中还提供更多的衍生格式。IEEE 754-1985 中规定的单精度浮点数和双精度浮点数分别称为 binary32 和 binary64。对每种格式都规定有最小的指数 emin 和最大的指数 emax。
可以用浮点数格式表达的有限值的范围取决于底数(b)、精度 (p) 和 emax。该标准一般将 emin 定义为1-emax：
 m 是 0 到 b^p-1 内的整数
例如，在 p=10 和 b=2 的情况下，m 就介于 0 到1023 之间。
 对给定的数，e 必须满足下式
1-emax<=e+p-1<=emax
例如，如果 p=24，emax=+127，则 e 的取值范围是-126 到 104。
这里需要清楚一点，即浮点表达式往往不是唯一的。
例如，0.02x10^1 和 2.00x10^(-1) 都代表相同的实数，即0.2。如果第一位数 d0 为 0，则称该数值被“标准化”。同时，对某个实数来说，可能不存在浮点表达式。例如，0.1 在十进制中是一个确定的值，但它的二进制表达式是一个小数点后 0011 的无穷循环。因此，0.1 无法用浮点格式确定地表达。表 1 给出了 IEEE 754-2008 规定的五种基本格式。


浮点计算的误差
因为要用固定位数来表达无穷数量的实数，因此在浮点计算中四舍五入是必不可少的。这就会不可避免地带来舍入误差，故需要有一种方法来衡量结果与采用无穷精度计算时的差距。我们来观察一下 b=10 和 p=4 的浮点格式。用这种格式，.0123456 可以表达成 1.234x10^(-2)。很明显这种表达方式在最后位置单位 (ulps) 发生了 .56 的差异。又如，如果浮点计算的结果是 4.567x10^(-2)，而采用无穷精度计算的结果是 4.567895，则最后差 .895 最后位置单位。
“最后位置单位”(ulps) 是规定这种计算误差的一种方法。相对误差是另一种用来衡量浮点数近似实数误差的方法。相对误差被定义为实数与浮点数之差除以实数的商。比如，将 4.567895 用浮点数表达为 4.567x10^(-2) 时的相对误差为 .000895/4.567895≈.00019。根据标准的要求，每个浮点数的正确计算结果的误差应不大于 0.5ulps。
浮点系统设计
在为数值应用开发设计的时候，很重要的一环是考虑输入的数据或者常数是否可以为实数。如果可以为实数，则在完成设计之前需要注意一些问题。需要检验来自数据集的数值在浮点表达式中是否很接近，有没有出现过大或者过小的情况。下面以二次方程求根为例。二次方程的根α 和 β 分别表达为下面两个等式：


很明显，如果 b^2 和 4ac 是浮点数，就会遇到舍入误差问题。如果 b^2 远大于 4ac，则 √(b^2 ) =|b|。如果b>0，β 将有解，但 α 的解将为 0，因为其分子中的各项相互抵消。另一方面，如果 b < 0，则 β 的解为 0，α 将有解。这是一种错误的设计，因为在浮点计算的作用下，根的值为 0，如果用于后续的设计，将很容易导致错误的输出。
防止这种错误抵消的方法之一是将 α 和 β 的分子和分母同时乘以 -b-√(b^2-4ac)，这样得到：


现在，在 b>0 的情况下，应使用 α’ 和 β 求出两个根的值。在 b<0 的情况下，应使用 α 和 β’ 求出两个根的值。在设计中很难察觉到这样的误差，原因在于对于整数数据集来说，常规公式会得出正确的结果，但在使用浮点数据的时候，常规公式会让设计返回错误的结果。
另一个有趣的实例是计算三角形面积。根据 Heron 公式，在已知三角形边长的情况下，三角形的面积 A 可以用下面的公式求出：


这里 a，b 和 c 是三角形的边长，s= (a+b+c)/2。
对一个一边的长度约等于另外两边长度之和的平三角形而言，有 a≈b+c，随即有 s≈a。由于两个相邻的数值会在 (s-a) 中相减，所以 (s-a) 可能会等于 0。或者，如果 s或者 a 出现舍入误差，就会导致显著的 ulps 误差。因此，最好根据数学家和计算机科学家 William Kahan (2)在 1986 年率先提出的建议，将此公式改写为下列格式：


我们看看分别用 Heron 的公式和 Kahan 的公式对相同的数值进行计算的结果如何。假定 a=11.0，b=9.0，c=2.04，则正确的 s 值为 11.02，正确的 A 值为1.9994918。但是，在采用 Heron 公式的情况下，计算得出的 s 值为 11.05，只有 3ulps 的误差，A 的值为3.1945189，出现巨大的 ulps 误差。在采用 Kahan 的公式的情况下，得出的 A 值为 1.9994918，小数点后七位完全与正确的结果吻合。
处理极小和极大的浮点数
极小的浮点数在相乘后会得出更小的结果。对设计人员来说，这可能会导致系统的下溢问题，最终的结果可能出错。这个问题在概率处理领域中尤为突出，比如计算语言学和机器学习。任何事件的概率一般都在 0 和 1 之间。能够得到完美整数结果的公式在处理非常小的概率时往往会导致下溢问题。避免下溢的方法之一是采用指数和对数。
在对很小的浮点数做相乘运算时，采用下面的公式非常有效：


这里指数和对数的底数必须相同。不过需要注意的是，如果a 和 b 都是概率，值都在 0和 1 之间，那么它们的对数会在-∞ 和 0 之间。如果只希望处理正概率，则将概率的对数乘以 -1，就可以正确地解读结果。采用对数值的另一项优势在于速度，因为相加比相乘的操作速度更快。
工程人员往往会遇到需要对非常大的连续数进行相乘的情况，比如：


K 的值在某些情况下非常容易溢出。处理这种情况的方法之一是用下式计算 k：


当然，这样计算相对较为复杂一点，但如果想要保证设计在计算这种极端状况时不会出错，这点代价还是值得的。
一些常见的浮点微妙之处
浮点算术并不严格遵循整数的数学规律。加法和乘法结合律对浮点数无效，减法和除法也是如此。造成这种情况的原因有两个。首先在顺序变动时会出现下溢或者上溢，其次是舍入误差问题。请看下面关于浮点数 a，b 和 c的几个例子。


右边的表达式在 a 和 b 相乘时可能出现上溢或者下溢。另外，由于 a 和 b 的乘积先计算，它的舍入误差与左边 b 和c 乘积的舍入误差不同，结果导致两个表达式得出的结果不同。


左边的 (b-c) 表达式中，如果两个数的浮点表示足够接近，会使该项结果为 0。在这种情况下 a=c。(b-c) 的舍入误差会导致结果与 b 不同。


由于舍入误差，表达式 (1.0+c) 的结果会导致值发生改变，最终与第一个表达式的结果不同。


同理，c/10.0 的舍入误差与 cx0.1 的不同，导致 b 的值不同。另外，数值 10.0 可以用 2 为底数表达为确切的表达式，但数值 0.1 用 2 为底数不能表达为确切的表达式。数值0.1 在表达为二进制表达式时，会出现 0011 的无穷循环。这种不能用特定的底数进行有限表达的数被称为该底数的非确定数，而能够进行有限表达的，则被称为确定数。
假定 c 可以一直被表达为确定数，如果 c 在一种情况下与确定数相除，在另一种情况下与非确定数相乘，则得到的结果是不一样的。而且由于 c 本身也可以是确定数或不确定数，取决于之前的计算方法，所以两个表达式之间的等价性无法成立。


这个等式在 x 为 0、无穷或者非数值的情况下是不成立的。在执行无效数学运算的情况下会出现 NaN，比如开负数的平方根或者用 0 除以 0。在出现 NaN 的情况下，程序员负责决定下一步怎么做。应注意根据 IEEE 754-2008 的规定，NaN 可以用指数 emax+1 和非零有效数表达为浮点数。因为有效数中可以放入由系统决定的信息，故结果将是一系列NaN 而非唯一的 NaN。因此在处理浮点代码时，请勿用 1.0替换 x/x，除非能够确认 x 不是零、无穷或者 NaN。


在 x 为无穷或者 NaN 的情况下该等式不成立。因此在不能确定 x 取值的情况下请勿用 0.0 替换 x-x。工程人员往往在编译器中采用这种替换技巧作为优化策略。但是正如所看到的，这样的优化可能导致错误的结果。


上面的等式在 x 为 -0 的时候不成立（注意可以有 +0 和-0）。这是因为根据 IEEE 754-2008 规定的默认舍入规则，(-0) + (+0) =+0，经舍入后为偶数。
IEEE 标准定义 +0 和 -0 的目的是为了得出答案的符号。比如 4 x (+0) =+0，(+0)/(-4)=-0。为比较目的，可以令+0=-0，这样像 if (x=0) 这样的声明就不用关心 +0 和 -0 之间的差异了。
更深入的理解
除了 Kahan 的著作，想要更加深入地领会使用浮点数据的各种微妙之处，工程人员还可以阅读 Goldberg 博士的《What Every Computer Scientist Should Know AboutFloating-Point Arithmetic》A 版（加利福尼亚州 Mountain View，1992 年 6 月）。另外有大量优秀参考文件对 IEEE 754-2008 标准以及其他影响浮点设计的标准（比如 ISO/ IEC 9899:201x）进行了详细解读。(3)
像航空电子和安全嵌入式系统这样的电子行业领域要求开发高度精确且无误的硬件和软件。由于部分硬件需要处理大量的浮点数据，因此设计人员有必要掌握浮点算术与整数算术之间的差异。
参考资料
1. IEEE 754-2008 标准，IEEE 浮点算术标准
2. W.Kahan，《计算针状三角形的面积和角度》，未出版，草稿见http://www.cs.berkeley.edu/~wkahan/
3. N1548， ISO/IEC 9899:201x，编程语言——C