[译文]任意精度运算的使用原因及方法
2013-11-17 12:17:02

按: 最近需要使用任意精度计算, 查阅资料时发现一篇文章(原文链接)写得浅显易懂, 并有示例代码. 遂译之, 为来者鉴.

There's a credibility gap: We don't know how much of thecomputer's answers to believe. Novice computer users solve this problem byimplicitly trusting in the computer as an infallible authority; they tend tobelieve that all digits of a printed answer are signifcant. Disillusionedcomputer users have just the opposite approach; they are constantly afraid thattheir answers are almost meaningless.
                       --Donald E.Knuth

这里存在信任度差异: 我们不知道计算机给出的答案有多少是可信的. 计算机新手完全信任计算机, 把它当作绝不出错的权威, 他们相信打印出的答案中所有的数字都有意义; 大失所望的计算机用户恰恰相反, 他们始终担心自己答案几乎毫无意义.
                       --Donald E.Knuth (高德纳)

任意精度运算的使用原因及方法
Kaveh R. Ghazi, Vincent Lefèvre Philippe, Théveny Paul Zimmermann

当今的浮点运算绝大多数采用双精度进行, 具有53比特的有效数字(significand或mantissa,参看附注名词解释). 然而, 有些应用需要使用扩展双精度(64比特或更多), 四精度(113比特)或更高精度. 在发表于2001年Astronomical Journal的一篇文章中, Toshio Fukushima写道:"在计算机日益强大的今天, 源于数值积分的误差仍是研究复杂动力系统--如太阳系和一些地外行星的长期稳定性--的主要障碍."他还通过一个实例指出, 对太阳系模拟而言, 使用双精度导致的舍入误差积累起来能够超过1弧度! 此外, 对运行于直升机或核反应堆电子控制系统的浮点程序进行静态分析, 使用任意精度也很必要.

假如我们要求出如下常数的10位有效数字(此例将贯穿本文始终),

$173746a + 94228b - 78487c$
其中$ a = sin(10^{22}), b = log(17.1), c =exp(0.42)$

在这个简单的例子里, 没有输入误差, 因为所有数字的值都精确知道, 具有无限精度.

下面是我们的第一个C语言程序,

#include <stdio.h>
#include <math.h>
int main (void)
{
   double a = sin (1e22), b = log(17.1);
   double c = exp (0.42);
   double d = 173746*a + 94228*b -78487*c;
   printf ("d = %.16en", d);
   return 0;
}

在64位Core2机器上, 操作系统Fedora10, 利用GCC 4.3.2(GNU libc 2.9)进行编译运行(下文测试平台相同), 我们得到

d = 2.9103830456733704e-11

这个结果完全错误, 因为正确结果为$-1.341818958 cdot 10^{-12}$. 我们可以将上面程序中的double换为long double以使用扩展双精度(64比特有效数字), 并将sin(1e22)换为sinl(1e22L),log换为logl, exp换为expl,%.16e换为%.16Le(注1),这样我们会得到

d = -1.3145040611561853e-12

这个新结果几乎错得和第一个一样离谱. 很显然我们计算时使用的精度不够.

注1: 在ARM计算机上,long double与双精度相同; 在Power PC上, long double相应于双双计算(参看附注); 在Solaris,long double对应于四精度.

1 哪里会出错?

在上面的例子中有些地方能够导致错误.

首先, 常数如1e22,17.1或0.42等可能无法以二进制精确表示. 对常数1e22,假定编译器能将它转换为正确的二进制数, 并符合IEEE 754的要求(参看附注IEEE 754), 它就能够在双精度下精确表示, 所以不存在这个问题. 然而,17.1并不能以二进制准确表示, 和它最接近的双精度值为$2406611050876109 cdot 2^{-47}$, 这与真值相差$1.4 cdot 10^{-15}$. 对0.42也存在类似的问题.

其次, 对一个数学函数(如sin)和一个浮点输入(如$x=10^{22}$),真值$sin x$通常并不能以双精度准确表示. 我们最多能将其舍入为最接近的双精度数$y=-7675942858912663 cdot 2^{-53}$, 误差$y-sin x$约为$6.8cdot 10^{-18}$.

再次, IEEE 754不仅没有规定对数学函数如$sin,log, exp$要进行正确的舍入, 就连精度也没有要求, 这导致计算结果完全取决于所用的平台[3].然而, 虽然IEEE1985版本没有提及这些数学函数, 2008版本中却将正确舍入列为推荐做法. 这样我们计算的$a,b, c$可能会和使用正确舍入方法得到的结果相差一些ulps(units inlast place, 参看附注名词解释).

在我们使用的测试平台上, 无论优化与否(参看第三节), 对由十进制输入舍入为相应二进制数的$x$,三个函数值都进行了正确的舍入.

最后, 在求$173746a+94228b-78487c$值的过程中出现抵消. 假定从左到右进行计算, $ 173746a+94228b $和$ 78487c $分别舍入为
$begin{split}
x&=1026103735669971·2^{-33}≈119454.19661583972629 \
y&=4104414942679883·2^{-35}≈119454.19661583969719
end{split}$
根据Sterbenz定理(参看附注名词解释), 计算$x-y$不存在舍入误差. 但最终结果的精确度明显地由计算$x$和$y$的舍入误差决定, 约为$ 2^{-36}≈1.5·10^{-11} $.由于要求的真值与此误差大小相近, 这样就导致了我们的最终计算结果d完全错误.

2 GNU MPFR库

论及任意精度, 我们指的是用户能选择每个运算的计算精度(也被称为多精度, 因为这意味着大的有效数字被分储在几个机器字中, 现代计算机一个字中最多能储存64比特, 约20位数字). 有一些程序和库能够用于任意精度计算, 大多数计算机代数系统, 像Mathematica,Maple和sage都支持任意精度计算. 本文中我们将重点介绍GNUMFPR, 一个用于任意精度浮点运算的C语言库(其他语言请参看附注其他语言).MFFR的特点在于运算时能够保证正确的舍入. 我们前面的程序可利用MPFR改写为:

#include <stdio.h>
#include <stdlib.h>
#include "mpfr.h"
int main (int argc, char *argv[])
{
   mp_prec_t p = atoi (argv[1]);
   mpfr_t a, b, c, d;
   mpfr_inits2 (p, a, b, c, d,(mpfr_ptr) 0);
   mpfr_set_str (a, "1e22",10, GMP_RNDN);
   mpfr_sin (a, a, GMP_RNDN);
   mpfr_mul_ui (a, a, 173746, GMP_RNDN);
   mpfr_set_str (b, "17.1",10, GMP_RNDN);
   mpfr_log (b, b, GMP_RNDN);
   mpfr_mul_ui (b, b, 94228, GMP_RNDN);
   mpfr_set_str (c, "0.42",10, GMP_RNDN);
   mpfr_exp (c, c, GMP_RNDN);
   mpfr_mul_si (c, c, -78487, GMP_RNDN);
   mpfr_add (d, a, b, GMP_RNDN);
   mpfr_add (d, d, c, GMP_RNDN);
   mpfr_printf ("d =%1.16Ren", d);
   mpfr_clears (a, b, c, d, NULL);
   return 0;
}

这个程序的输入参数为计算精度. 取$p=53$,运行结果为

d = 2.9103830456733704e-11

这和我们使用双精度得到的结果完全一样.$p=64$时, 结果为

d = -1.3145040611561853e-12

与扩展双精度结果一致.$p=113$时, 相应于IEEE-754的四精度, 结果为

d = -1.3418189578296195e-12

此结果与真值完全一致.

3 常数合并

在一个给定的程序中, 若一个表达式是常数, 如$3+ (17 times 42)$, 很可能在编译时就被替换为它的计算值. 这对浮点值也是一样, 却另有难处: 编译器需要能够确定使用的舍入模式. 这种编译器所做的替换被称为常数合并[4].使用正确舍入的常数合并, 所产生的常数只取决于目标平台的浮点格式, 而与创建平台的处理器, 系统, 数学库等无关. 这既保证了正确性(产生的常数具有正确的精度和舍入模式)也保证了可重复性(去除了平台相关性).GCC 4.3版本使用MPFR对内部数学函数, 如sin,cos, log, sqrt, 实现常数合并. 考虑下面的例子

#include <stdio.h>
#include <math.h>
int main (void)
{
   double x = 2.522464e-1;
   printf ("sin(x) = %.16en",sin (x));
   return 0;
}

使用GCC 4.3.2, 如果编译时关闭优化(使用-O0),我们得到结果2.4957989804940914e-01. 若打开优化(-O1),我们将得到2.4957989804940911e-01. 结果为什么不同? 使用-O0时, 表达式sin(x)由数学库计算(使用的是GNUC语言库, 也称为GNUlibc或glibc). 使用-O1时,GCC将表达式sin(x)识别为常数, 利用最接近舍入模式, 调用MPFR计算sin(x)的值, 并直接将其替换为正确的舍入值(注2).正确的值是使用-O1时的值. 需要注意的是, 即便GNUC库并没有对上面的例子进行正确地舍入, 它对大多数值都能进行正确地舍入(在没有扩展精度的计算机上), 正如IEEE 754-2008所推荐的. 在将来, 我们希望对每个输入和函数都能进行正确的舍入.

注意: 在x86处理器上, 若GNUC语言库使用x87处理器的fsin实现, 对$x= 0.2522464$能够给出正确的舍入值. 然而, 这只是巧合, 因为在0.25及更大的$10^7$个双精度数字中, 有2452个fsin不能给出正确的舍入.

注2: 当使用-O1编译时, 我们甚至可以忽略链接数学库. 运行gcc-O1 test.c根本就没有使用数学库. 与此相对, 运行gcc-O0 test.c将会得到一个编译错误, 而gcc-O0 test.c -lm则正常. 这说明使用-O0时需要数学库. 取消常数合并和其他对内部函数的优化, 可使用gcc-fno-builtin, 或使用更具体的gcc -fno-builtin-sin取消对sin函数的优化.

4 结论

在本文中, 我们可以看到, 使用具有53比特有效数字的双精度变量能够导致远低于53比特精度的最终结果. 导致精度丢失的可能原因包括舍入误差, 数字抵消, 二进制-十进制转换误差, 低数值品质的数学函数,...使用任意精度运算, 例如GNUMPFR库, 能够提高最终计算结果的精确度. 更重要的是,MPFR能够对数学函数进行正确地舍入, 这有助于提高在不同处理器, 编译器和操作系统下进行的浮点运算的可重复性, 正如我们在上面GCC常数合并的例子中所展示那样.

参考文献
[1] Fousse, L., Hanrot, G., Lefèvre, V., Pélissier, P., and Zimmermann, P.MPFR: A multiple-precision binary floating-point library with correct rounding.ACM Trans. Math. Softw. 33, 2 (2007), article 13.
[2] IEEE standard for floating-point arithmetic. ANSI-IEEE standard 754-2008,2008. Revision of ANSI-IEEE Standard 754-1985, approved June 12, 2008: IEEEStandards Board.
[3] Lefèvre, V. Test of mathematical functions of the standard C library.http://www.vinc17.org/research/testlibm/.
[4] Wikipedia. Constant folding. http://en.wikipedia.org/wiki/Constant_folding.

附注

IEEE 754标准
IEEE 754是一个广泛使用的关于浮点数表示和操作的标准(你的计算机每天都在使用它, 你甚至都没注意到). 它非常重要, 因为其中定义了浮点数格式--能够使两台计算机互相交换浮点值而不损失任何精度--并要求对算术运算正确舍入, 这保证了同样的程序在不同计算机上能够得到完全相同的结果(注3).IEEE 754在1985年首先通过, 并在2008年进行了修订[2].我们将修订后的这个版本称为IEEE 754-2008. 它定义了基本格式(用于运算)和交换格式(用于不同实现间的数据交换). 共有五种基本格式:binary32, binary64, binary128, decimal64和decimal128.binary32和binary64分别促生了单精度和双精度, 通常对应于ISOC语言的float和double数据类型. 十进制格式是IEEE754-2008中新定义的, GCC只提供了初步的支持. 例如,decimal64在GCC中被标记为_Decimal64,与目前TR 24732草稿中关于C语言十进制浮点算术的规定一致(注4).把我们的示例程序改写成:

#include <stdio.h>
#include <math.h>
int main (void)
{
   _Decimal64 a = sin (1e22);
   _Decimal64 b = log (17.1);
   _Decimal64 c = exp (0.42);
   _Decimal64 d =173746*a+94228*b-78487*c;
   printf ("d = %.16en",(double) d);
   return 0;
}

得到结果

d = 0.0000000000000000e+00

(注意必须把最终结果d转换为double,因为GNU C语言库还不支持十进制格式输出).

IEEE 754-2008要求对四类基本算术运算(+,-,×,÷), 平方根和基数转换(如将十进制字符串读入为二进制格式, 或将二进制格式输出为十进制字符串)进行正确地舍入. 这意味着计算结果好像应该以无限精度进行, 然后再按相应的舍入模式进行舍入.IEEE 754-2008明确了五种舍入模式(或属性): roundTowardPositive, roundTowardNegative, roundTowardZero, roundTiesToEven和roundTiesToAway.

注3: 需要一些我们在这里忽略的条件
注4:http://www.open-std.org/jtc1/sc22/wg14/www/projects

扩展精度和Linux
传统的32位x86处理器的浮点运算单元能够设定为对计算结果以双精度(53比特有效数字)或扩展精度(64比特有效数字)进行舍入. 大多数操作系统, 如FreeBSD,NetBSD和Microsoft Windows, 默认都选择以双精度进行舍入. 而Linux默认以扩展精度进行舍入. 这是一个坏的选择, 详细的原因请参看http://www.vinc17.org/research/extended.en.html.

名词解释

基数, 有效数字和指数
若$x$是一个浮点数, 精度为$p$,基数为$beta$, 它能被写为$x= pm 0.d_1d_2 . . . d_p cdot beta^e$, $s = pm 1$是$x$的符号.$m = 0.d_1d_2 . . . d_p$是$x$的有效数字,$e$是$x$的指数. 注意, 这种表示方法并不唯一, 除非要求$d_1$为非零值. 此外, 有效数字可以采用不同的约定, 并因此而改变指数的值. 例如,IEEE 754-2008使用$m = d_1.d_2 . . . d_p$约定, 使得其指数比我们定义的小1. 它也使用了第三个约定,$m$是一个整数.

Unit in lastplace
若$x = pm 0.d_1 d_2 . . . d_p cdotbeta^e$是一个浮点数, 定义$mathrm{ulp}(x)$为$x$最小有效数字的权重, 即$beta^{e-p}$.(注意$mathrm{ulp}(x)$的值并不取决于$(s,m,e)$约定)

Sterbenz定理
若$x$和$y$是两个浮点数, 精度为$p$,并满足$y/2 ll x ll 2y$, 则$x-y$能够以精度$p$准确表示. 作为推论, 计算$x-y$时没有舍入误差.

双双算术
将一个实数$r$近似为两个双精度数的和,$x+y$. 若$x$是$r$的最近舍入,$y$是$r-x$的最近舍入, 双双算术的精度是使用单个双精度书的两倍.

其他语言
除C和C++外, 其他一些语言也支持任意精度浮点运算.Perl, Python, Haskell, Lisp和Ursala都有MPFR的接口(详情请参看mpfr.org).在Fortran中使用MPFR并不容易, 因为必须首先将MPFR的数据类型转换为Fortran支持的. 但是, 如果你需要计算的只是双精度数的指数函数的正确舍入, 很容易, 参看http://www.loria.fr/~zimmerma/mpfr/fortran.html.

如何在Fortran中使用MPFR
Credit: thanks to Olivier Cessenat (CEA) for his help.

首先, 创建一个C文件,crexp.c

#include <mpfr.h>
double exp_(double *x) {
 double y ;
 mpfr_t z ;
 mpfr_init2 (z, 53);
 mpfr_set_d(z, *x, GMP_RNDN) ;
 mpfr_exp(z, z, GMP_RNDN) ;
 y = mpfr_get_d(z, GMP_RNDN) ;
 mpfr_clear(z) ;
 return(y) ;
}

然后, 创建一个Fortran程序,toto.f90

program toto
 real(kind=8) :: x, y
 real(kind=8),external :: exp
 x = 1.D0
 y = exp(x)
 print*,'y=',y
end program toto

编译

gcc -c crexp.c -o crexp.o
gfortran -c toto.f90 -o toto.o
gfortran -o a.out toto.o crexp.o -lmpfr

运行

./a.out

译后记:

查看MPFR, 发现awk4.1版本已经支持MPFR, 而且使用非常简单. 这对我这个awk的深度用户来说, 值得一记, 也让我对awk的喜爱更添了几分.