2013-11-17 12:17:02
按: 最近需要使用任意精度计算, 查阅资料时发现一篇文章写得浅显易懂, 并有示例代码. 遂译之, 为来者鉴.
There’s a credibility gap: We don’t know how much of the computer’s answers to believe. Novice computer users solve this problem by implicitly trusting in the computer as an infallible authority; they tend to believe that all digits of a printed answer are signifcant. Disillusioned computer users have just the opposite approach; they are constantly afraid that their 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语言程序,
# Language: 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 = %.16e\n", d);
return 0;
}
在64位Core2机器上, 操作系统Fedora 10, 利用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.8 \cdot 10^{-18}$.
再次, IEEE 754不仅没有规定对数学函数如 $\sin, \log, \exp$ 要进行正确的舍入, 就连精度也没有要求, 这导致计算结果完全取决于所用的平台[3]. 然而, 虽然IEEE 1985版本没有提及这些数学函数, 2008版本中却将正确舍入列为推荐做法. 这样我们计算的 $a, b, c$ 可能会和使用正确舍入方法得到的结果相差一些ulps(units in last place, 参看附注名词解释).
在我们使用的测试平台上, 无论优化与否(参看第三节), 对由十进制输入舍入为相应二进制数的 $x$, 三个函数值都进行了正确的舍入.
最后, 在求 $173746a + 94228b -78487c$ 值的过程中出现抵消. 假定从左到右进行计算, 和 $173746a +94228b$ 舍入为 $x = 1026103735669971 \cdot 2^{-33} ≈ 119454.19661583972629$, $78487c$ 舍入为 $y = 4104414942679883 • 2^{-35} =119454.19661583969719$. 根据Sterbenz定理(参看附注名词解释), 计算 $x-y$ 不存在舍入误差. 但最终结果的精确度明显地由计算 $x$ 和 $y$ 的舍入误差决定, 约为 $2^{-36} ≈ 1.5•10^{-11}$. 由于要求的真值与此误差大小相近, 这样就导致了我们的最终计算结果d完全错误.
2. GNU MPFR库
论及任意精度, 我们指的是用户能选择每个运算的计算精度(也被称为多精度, 因为这意味着大的有效数字被分储在几个机器字中, 现代计算机一个字中最多能储存64比特, 约20位数字). 有一些程序和库能够用于任意精度计算, 大多数计算机代数系统, 像Mathematica, Maple和sage都支持任意精度计算. 本文中我们将重点介绍GNU MFPR, 一个用于任意精度浮点运算的C语言库(其他语言 请参看附注其他语言). MFFR的特点在于运算时能够保证正确的舍入. 我们前面的程序可利用MPFR改写为:
# Language: c
#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.16Re\n", 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, 实现常数合并. 考虑下面的例子
# Language: c
#include < stdio.h>
#include < math.h>
int main (void)
{
double x = 2.522464e-1;
printf ("sin(x) = %.16e\n", sin (x));
return 0;
}
使用GCC 4.3.2, 如果编译时关闭优化(使用-O0), 我们得到结果2.4957989804940914e-01. 若打开优化(-O1), 我们将得到2.4957989804940911e-01. 结果为什么不同? 使用-O0时, 表达式sin(x)由数学库计算(使用的是GNU C语言库, 也称为GNU libc或glibc). 使用-O1时, GCC将表达式sin(x)识别为常数, 利用最接近舍入模式, 调用MPFR计算sin(x)的值, 并直接将其替换为正确的舍入值(注2). 正确的值是使用-O1时的值. 需要注意的是, 即便GNU C库并没有对上面的例子进行正确地舍入, 它对大多数值都能进行正确地舍入(在没有扩展精度的计算机上), 正如IEEE 754-2008所推荐的. 在将来, 我们希望对每个输入和函数都能进行正确的舍入.
注意: 在x86处理器上, 若GNU C语言库使用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比特精度的最终结果. 导致精度丢失的可能原因包括舍入误差, 数字抵消, 二进制-十进制转换误差, 低数值品质的数学函数, …使用任意精度运算, 例如GNU MPFR库, 能够提高最终计算结果的精确度. 更重要的是, 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: IEEE Standards 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分别促生了单精度和双精度, 通常对应于ISO C语言的float和double数据类型. 十进制格式是IEEE 754-2008中新定义的, GCC只提供了初步的支持. 例如, decimal64在GCC中被标记为_Decimal64, 与目前TR 24732草稿中关于C语言十进制浮点算术的规定一致(注4). 把我们的示例程序改写成:
# Language: c
#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 = %.16e\n", (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 last place
若 $x = \pm 0.d_1 d_2 . . . d_p \cdot \beta^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
# Language: 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
# Language: fortran
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, 发现awk 4.1版本已经支持MPFR, 而且使用非常简单. 这对我这个awk的深度用户来说, 值得一记, 也让我对awk的喜爱更添了一层.
