科学计算当中的高精度浮点运算

1、科学计算当中的高精度浮点运算David H. Bailey 2005年1月25日 摘要：目前，IEEE的64位浮点运算在大多数科学应用当中足够精确。然而，对于 快速增长的重要科学计算应用中需要更高水平的数值精度其中的一些应用 需要大约两倍于现在的水平其他的需要四倍;而还有一些人为了获得有意义 的数值结果需要数百或更多位数新的高精度的软件包包括用来减少转换工 作的高级语言翻译模块为这样的计算提供了便利在不同领域里包括量子理 论，气候建模与数学实验，这些活动已经产生了一些有趣的新的科学成骐 中一些将在本文中所描述这些进展表明，在未来，用于科学计算的数值精度 同程序设计的算法和数据结构一样重要。关键

2、词：高精度算法，数学实验气候模拟，量子理论，计算化学，计算物理 说明：这项工作由计算技术研究，离散数学数学，信息，美国计算科学和能源部 支持，合同号 :DE-AC03-76SF000981、序言几乎所有的现代计算机系统，从个人计算机到最大的超级计算机，实现了 64位浮点运算标准，提供3个尾数位，或大约6位十进制精度。对于大多数科 学应用来说，这是绰绰有余的，以及一些应用程序，如常规试验数据的处理， 即使是32位的标准彳主彳主也恩能提供足够的精度。然而，对于一个迅速扩大的应用机构64位的算术将不再是足够的。这些 范围从一些有趣的新的数学计算到在高性能并行超级计算机上进行的大型物理 模拟。在这些应

3、用中，部分包括敏感的计算数值的常规代码，使用常规算法产 生结果的准确性有问题。这些不准确地结果可能反过来导致其他错误，如在一 个条件分支上错误的路径。一些困难，这种类型的可能是通过改变基本的算法或改变某些操作的顺序 从而被修改，其中一些操作可移执行。但在其他情况下，这种修正是有效的， 其中一些问题，是正常水平的两倍，大约-digit精度是足够的。在其他情况 下，四倍的正常水平是必要的或大约六十四位数的精度在其他一些情况下， 包括从物理和数学的应用，更高水平的精度是必要的，从数百到数千位。重要的是要记住这里非常少数的目前参与技术计算的科学家和工程师谁是 有严格的数值分析背景中。更重要的是，这种稀

4、缺的数值技术有可能恶化，如 果没有改善的话，在未来，部分原因是由于缺乏对数值计算有兴趣的学生和年 轻的研究人员22。因此，虽然有些人可能认为，敏感的数值计算可以通过 使用不同的算法或编码技术，而在实践中这种变化是非常容易出错而且人力成 本非常昂贵。换句话说，使用高精度算法克服这些困难，通常将更容易，更便 宜和更可靠，即使其他的补救措施在理论上是可行的。本文将着重于高精度的浮点数而不是整数运算。然而，应该注意的是，高 精度（即多字）整数运算在所有的舞台上本身就是一个非常有趣的，有众多的 应用在纯数学和应用数学。例如，当登录到一个安全网站购买书或电脑配件， 浏览器软件进行高精度的整数计算交流信用卡

5、号码和其他安全信息。与此密切 相关的技术是在大整数分解算法进行的研究，其中最近我们已经取得了几个显 著的进展。一个很好的参考这是最近出版的兰德尔和波默朗斯的一本书。在crandall-pomeranc书中描述的jnathe-matic的算法的实现可从 /software.h找找 到。2、高精度软件进行高精度浮点运算的软件包在早期的计算当中就已经实现。例如，李察布兰特写的一个软件包0世纪70年代以来就开始被使用14。然而，许多这些软件 包需要重写。很难作出这样的改变，以及难以调试产生的代码，这阻止了只有 少数科学家能够利用这种软件。在过去几年中，高精度的软件包已产生，包括 高级语言的接口，使这种

6、转换相对不那么困难。这些软件包通常利用自定义的 数据类型，运算符重载的特点，现在在语言如I+和Fortran-90E也可以，转换 很方便。更先进的高精度计算在商业软mathematical Maple中是可以的这些产 品以非常自然的方式将任意精度算术，在纯数学和应用数学中，延伸到了一个 范围广泛的先进功能。这些软件包在一些方面中是非常有用的，但他们一般没 有专业高精度的计算软件包快，同时rtran-9或C+编程中提到，他们没有 提供一种手段，利用高精度的计算设施来转换现有的科研项目。一些高精度的计算软件包例子，在互联网上免费提供的有下面这些（按字 母顺序列出）ARPREC, DDFUN, QD

7、 and MPFUN在作者的网站上可以得到： HYPERLINK /dhbailey/mpdist /dhbailey/mpdist or from HYPERLINK . ARPREC .这个软件包包括例程执行任意高虐但需设定）的精度水平的运算， 包括许多代数和超越函数高级语言的接口对于+和Fortran-9是可行的，配 套整数和复杂数据类型。该软件Xiaoye S. Li and Brandon Thompson.编写。 DDFUN这个软件包在所有ortran-9（环境中，提供了双倍（约31位数）算法。高级语言的接口支持实数，整数和复杂数据类型。这个软胶包的速度远远超过 正在使用的任意精度

8、或在应用中的个数字足够够的iuad-doubl软件，而且彳主 彳主比在一些语言系统中使用粉al*1（数据类型更快。FMLIB和FMZM90. FMLI是低层次的多精度库，并提供高级语言接口。同时 FMZM90为实数，整数和复杂数据类型提供高FSrtran-90言接口。这个软件 是由戴维斯密的忠臣大学编写。在以下网站可以得到： HYPERLINK /dmsmith/fmlib.h /dmsmith/fmlib.h k 0有74个四字节的词。这里面只有两个能够在原始条件下发生循环环费曼图 中：V (-1) j cos(2 兀 k / 3)V := Req (O2w w )= 乙 j k 0j *从

9、图表中我们发现在较少循环的费曼图中，剩余项的常数结果也减少。这十种情况如图3所示。在这些图中，点粒子的静止质量显示为零。Q,该公已被发现，相应的常数在图在这里，常数C = sin(兀 k /3) / k 2在每一种情况下，这些公式后来被证明是正确的。但直到“实验”的计算 结果被发现时，没有人有任何坚实的理由相信这种关系可能存在。实上，“实 验”这个形容词是比较合适的，因为PSLQ计算与常规实验室实验发挥了作 用完全类似（只是一个至关重要）图3：十个四面体的例子13、结论我们在这里提出一个在现代科学计算中迅速扩大使用高精度算法简短的调 查。值得注意的是，所有这些例子已经出现在过去的十年。因此，我

10、们可以看 到科学计算这个新时代的诞生，在计算精度要求的计算是中程序设计的算法和 数据结构是非常重要的。VI| K(3)+3C(4)巧且=6 一51弓N=6-号-8U=6一火W=6=3) 京(4) 一=6(3)一射(4) 63V4A=6-)-63=6(3) 14(4) 16*=6=3) 舞0(4) + 沁2 一 16/=6(3) - 13(4) - 8U - 4C2图4：由PSLQ公式发现的图的十个例子在一个方面，这也许不是巧合，在高精度的计算已经出现在同一时期，一 个共同的努力在许多科学计算的高度并行和分布式系统上已实施。这些系统比 以前提供了更大规模的运行，极大地放大的数值困难。此外，即使在

11、单处理器 系统中，内存和性能已经大大扩展（在半导体技术中穆尔定律的礼物），使更 大的计算与几年前相比是可行的。这很可能是在这些系统正在执行的其他计算显著的数值困难，但那些使用 这些代码的没有意识到这些困难的程度。随着使用高精度算法相对简单的转换 程序软件的出现，现在可以定期运行大型科学计算与双或更高精度的算法，以 确定如何保证结果是真实性。事实上，目前可用于高精度算法的软件工具可以用来在每个阶段的数值错 误管理：（1）定期检测科学的代码，看看它是否变变为了敏感数值;2）（确 定有多少位的结果（中间或最终）是可靠的；3）（用几行代码降低数值困难；（4）发现整流数值困难。本文描述的软件和其他地方在

12、短期来看将帮助科学家，从长远来看，至关 重要的是，计算机厂商意识到需要提供硬件和软件支持高精度算法。一个国际 委员会已经形成更新EEE-754浮点运算标准；128位的支持是标准草案的一 个关键条款。这个委员会也正在考虑建立标准的任意精度算术。这些都是受欢 迎和需要及时发展的。鸣谢作者承认许多同事在准备这篇文章的贡献，括:Ed Baron, BenBarrowes,Jonathan Borwein; Richard Crandall, Chris Ding, Alexei Frolov, WilliamKahan, Peter Nugent, Michael Strayer and Zong-C

13、hao Yan.参考文献：David H. Bailey,“Integer Relation Detectiorl,Computing in Science and Engineering Jan-Feb., 2000, pg. 24-28.David H. Bailey, Peter B. Borwein, and Simon Plouffe, “On the Rapid Computation of Various Polylogarithmic Constants, Mathematics of Computation, vol. 66, no. 218 (Apr 1997), pg.

14、903-913.David H. Bailey and David Broadhurst, “Parallel Integer Relation Detection: Techniques and Applications” Mathematics of Computation vol. 70, no. 236 (2000), pg. 1719-1736.David H. Bailey and David J. Broadhurst,“A Seventeenth-Order Polylogarithm Ladder, HYPERLINK /dhbailey/dhbpapers/ladder.p

15、df /dhbailey/dhbpapers/ladder.pdf (1999).David H. Bailey and Richard E. Crandal“,On the Random Character of Fundamental Constant Expansions Experimental Mathematics vol 10, no. 2 (June 2001), pg. 175-190.David H. Bailey and Richard E. Crandall, “Random Generators and Normal Numbers, Experimental Mat

16、hematics vol. 11, no. 4 (2004), pg. 527)46.D. H. Bailey and A. M. Frolov, “Universal Variational Expansion for High-Precision Bound-State Calculations in Three-Body Systems. Applications to Weakly-Bound, Adiabatic and Two-Shell Clouster Systems, Journal of Physics B, vol. 35, no. 20 (28 Oct 2002), p

17、g. 42870-4298.David H. Bailey, R. Krasny and R. Pelz, “Multiple Precision, Multiple Processor Vortex Sheet Roll-Up Computation, Proceedings of the Sixth SIAM Conference on Parallel Processing for Scientific Computing1993, SIAM, Philadelphia, pg. 5-56.David H. Bailey and Xiaoye S. Li,“A Comparison of

18、 Three High-Precision Quadrature Schemes, manuscript available at HYPERLINK /dhbailey/dhbpapers/quadrature.pdf /dhbailey/dhbpapers/quadrature.pdf (2004).David H. Bailey and Sinai Robins,“Highly Parallel, High-Precision Numerical Integration, manuscript available at HYPERLINK /dhbailey/dhbpapers/quad

19、parallel.pdf /dhbailey/dhbpapers/quadparallel.pdf (2004).Edward Baron and Peter Nugent, personal communication, Nov. 2004.B. E. Barrowes, K. O Neill, T. M. Grzegorczyk and J. A. Kong, “Asymptotic Expansions of the Prolate Angular Spheroidal Wave Function for ComplexeSParameter Studies in Applied Mat

20、hematics, to appear.Jonathan M. Borwein and David H. Bailey, Mathematics by Experiment: Plausible Reasoning in the 21st Century.Richard P. Brent, “A Fortran Multiple-Precision Arithmetic Package, ACM Transactions on Mathematical Softwre, vol. 4, no. 1 (Mar 1978), pg. 5-70.David J. Broadhurst, “Massi

21、ve 3-loop Feynman Diagrams Reducible to SC Primitives of Algebras of the Sixth Root of Unity, preprint, March 1998, to appear in European Physical Journal C. Available from HYPERLINK /abs/hepth/9803091 /abs/hepth/9803091.Russel E. Caflisch, “Singularity Formation for Complex Solutions of the 3D Inco

22、mpressible Euler Equations,Physica D vol. 67, no. 1 (1993), pg.-18.Richard Crandall and Carl Pomerance,Prime Numbers: A Computational Perspective, Springer-Verlag, New York, 2001.A. M. Frolov and D. H. Bailey, “Highly Accurate Evaluation of the Few-Body Auxiliary Functions and Four-Body Integrals Jo

23、urnal of Physics B, vol. 36, no. 9 (14 May 2003), pg. 1857-1867.Wayt Gibbs,“A Digital Slice of P, Scientific American vol. 288, no. 5 (May 2003), pg. 23-24.P. H. Hauschildt and E. Baron, “The Numerical Solution of the Expanding Stellar Atmosphere Problem, Journal Computational and Applied Mathematics vol. 109 (1999), pg. 41-63.Yun He and Chris Ding, “Using Accurate Arithmetics to Improve Numerical Repro-ducibility and Stability in Parallel Applications Journal of Supercomputing, vol. 18, no. 3 (Mar 2001), pg 259-277.William Kahan, personal communication, Nov. 2004.William Kahan and