（数学专业论文）多项式函数的高精度计算及舍入误差分析研究.pdf

国防科学技术大学研究牛院学何论文 摘要 在计算机技术高速发展的今天，多项式的函数值及其导数值评估在越来越多的领 域发挥着重要作用，几乎所有的函数在计算机运算过程中都是以多项式的形式来参与 运算。与此同时，针对大规模科学计算和病态问题中舍入误差的累积可能使算法最终 结果失真的问题，科研工作者需求高精度的可靠算法。因此，对计算多项式函数值和 导数值的高精度算法的需求就显得尤为重要和迫切。 本学位论文应用无误差变换技术针对幂指数基多项式一阶导数值并i c h e b y s h e v 基多 项式函数值的计算提出了新的高精度算法，取得了较为满意的结果，主要内容如下： 与计算多项式函数值相比，计算多项式导数值不但要考虑计算过程中产生的舍入 误差，还要考虑参数的累计误差。本文通过分析参数自身带入的累积误差与舍入误差 之间的关系，提出了补偿h o m e r 导数算法。该算法将传统算法的结果误差转化为三个以 舍入误差为系数的多项式的和，并将其补偿给原算法的数值结果。理论分析和数值实 验均表明改进算法具有高精度性和可靠性。 针对c h e b y s h e v 基多项式函数值计算，本文在c l e n s h a w 算法的基础上，提出传统算法 的结果误差为一个以舍入误差为系数的c h e b y s h e v 基多项式。将这一多项式的函数值补 偿给原算法数值结果，即构成本文的补偿c l e n s h a w 算法。理论分析和数值实验均表明改 进算法比传统算法优越。 关键词：补偿h o m e r 导数算法；补偿c l e n s h a w 算法；幂指数基多项式；c h e b y s h e v 基 多项式；无误差变换；舍入误差；i e e e 7 5 4 标准 第1 页 国防科学技术大学研究生院学位论文 a b s t r a c t w i t ht h er a p i dd e v e l o p m e n to fc o m p u t e rs c i e n c ea n dt e c h n o l o g y , p o l y n o m i a la n di t sd e r i v a t i v e s p l a yak e yr o l ei nm o r ea n dm o r ea p p l i e df i e l d s ，i nm o s to fn u m e r i c a lc o m p u t a t i o ni nm a c h i n et h e f u n c t i o ni sa p p r o x i m a t e db yt h ep o l y n o m i a l m e a n w h i l e ，s i n c et h ep r o p a g a t i o no fr o u n d o f fe r r o r s t h r o u g h o u tal a r g e s c a l eo ri l l c o n d i t i o n e dc a l c u l a t i o nm i g h tl e a dt ot h ec a t a s t r o p h i cn u m e r i c a lr e s u l t ， t h es c h o l a r sn e e da c c u r a c ya n dv a l i d a t e da l g o r i t h mi nn u m e r i c a lc o m p u t a t i o nf i e l d s c o n s e q u e n t l y , a c c u r a c ye v a l u a t i o no fap o l y n o m i a la n di t sd e r i v a t i v ei su r g e n t l yr e q u i r e d t h i sd i s s e r t a t i o nm a i n l yc o n s i d e r sc o n s t r u c t i n gt w oa c c u r a c ya n dv a l i d a t e da l g o r i t h m sb a s e d o ne r r o r - f r e et r a n s f o r m a t i o n ，n a m e df o rc o m p e n s a t e dh o m e rd e r i v a t i v ea l g o r i t h ma n dc o m p e n s a t e d c l e n s h a wa l g o r i t h mr e s p e c t i v e l y o n ea l g o r i t h me v a l u a t e st h ef i r s td e r i v a t i v eo fap o l y n o m i a lr e p r e s e n t e di np o w e rb a s i s ，t h eo t h e ri sf o re v a l u a t i o no fap o l y n o m i a li nc h e b y s h e vf o r m t oap o l y n o m i a li np o w e rb a s i s ，t h ei m p r o v e da l g o r i t h mi sb a s e do nt h eh o m e rd e r i v a t i v e m e t h o d 。u t i l i z e se r r o r - f r e et r a n s f o r m a t i o nt e c h n i q u et on o t er o u n d o f fe r r o r si nt h er e c u r s i o n t h e e r r o rg e n e r a t e db yt h eo r i g i n a la l g o r i t h mi sp r o v e dt ob ee x a c t l yt h es u mo ft h r e ep o l y n o m i a l s ，t h e c o e f f i c i e n t so fw h i c ha r et h er o u n d i n ge r r o r s f i n a l l y , t h ei m p r o v e da l g o r i t h mg e t sac o r r e c t e dr e s u l t b ya d d i n gt h i ss u mt ot h eo r i g i n a lr e s u l ti nt h ee n d t h e o r e t i c a la n a l y s i sa n dn u m e r i c a le x p e r i m e n t s i l l u s t r a t et h a tt h ei m p r o v e da l g o r i t h mc a np r o v i d el a s t b i ta c c u r a c y t oap o l y n o m i a li nc h e b y s h e vf o r m ，t h ei m p r o v e da l g o r i t h mp r o p o s e st h a tt h ee r r o ro ft h e c l e n s h a wa l g o r i t h m sn u m e r i c a lr e s u l ti se q u a lt oap o l y n o m i a li nc h e b y s h e vb a s i sw i t hf l o a t i n g p o i n tc o e f f i c i e n t s ，a n dt h e s ec o e f f i c i e n t sa r eg e n e r a t e db yt h ec l e n s h a wa l g o r i t h m t h ei m p r o v e d a l g o r i t h mg e t sac o r r e c t e dr e s u l tb ya d d i n gt h i sn e wp o l y n o m i a lt ot h eo r i g i n a lr e s u l t t h e o r e t i c a l a n a l y s i sa n dn u m e r i c a le x p e r i m e n t si l l u s t r a t et h a tt h ei m p r o v e da l g o r i t h mi sa c c u r a c ya n dv a l i d a t e d k e yw o r d s ：c o m p e n s a t e dh o m e rd e r i v a t i v ea l g o r i t h m ，c o m p e n s a t e dc l e n s h a wa l g o - r i t h m , m o n o m i a lp o l y n o m i a l ，c h e b y s h e vp o l y n o m i a l ，e r r o r - f r e et r a n s f o r m a t i o n ，r o u n d - o f f e r r o r , i e e e - 7 5 4s t a n d a r d 第l i 页 国防科学技术大学研究牛院学位论文 第i i 页 独创性声明 本人声明所呈交的学位论文是我本人在导师指导下进行的研究工作及 取得的研究成果。尽我所知，除文中特别加以标注和致谢的地方外，论文中 不包含其他人已经发表和撰写过的研究成果，也不包含为获得国防科学技 术大学或其他教育机构的学位或证书而使用过的材料。与我一同工作的同 志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意。 学位论文题目：多项式函数的高精度计算及舍入误差分析研究 学位论文作者签名： 壬 拉日期：夕哆年，月日 学位论文版权使用授权书 本人完全了解国防科学技术大学有关保留、使用学位论文的规定。本人 授权国防科学技术大学可以保留并向国家有关部门或机构送交论文的复印 件和电子文档，允许论文被查阅和借阅；可以将学位论文的全部或部分内 容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、 汇编学位论文。 ( 保密学位论文在解密后适用本授权书。) 学位论文题目：多项式函数的高精度计算及舍入误差分析研究 学位论文作者签名：三壬宰拯日期：、j r ， - zl = 一一 作者指导老师签名：邀n 盈日期： 刁年，月日 神年| | 疑t sb ， 围防科学技术大学研究生院学位论文 第一章绪论 1 1 本文研究背景及现况 随着科学技术的进步，在越来越多的大规模科学计算中科学家们广泛的需要可靠 的高精度数值计算方法。目前应用较为广泛的多精度算法库有固定字长扩展如双倍双 精度( d o u b l e d o u b l e ) 和四倍双精度算法( q u a d d o u b l e ) ，用来模拟两倍和四倍的i e e e 7 5 4 双 精度，见参考文献 1 5 1 ；作为这类固定字长的算法推广和发展x b l a s 提出了扩展的混合 精度算法，见 5 2 矛n 2 6 1 ；另外，还有法国i n r i a 实验室n r e v o l 等人提出的m p f r ，m p f i ，见 文献 3 2 1 ，还有j v i g n e s 4 9 提出的c e s t a c ，d h b a i l e y 2 提出的a r p r e c ，以及r m o o r e 提 出的区间算法，见文献 3 1 1 和【3 4 】。此外，本论文中将要引用的补偿算法也是一类应用 较为广泛的有效算法。算法只使用i e e e 7 5 4 双精度来实现高精度和及内积的运算，早 期研究学者有t j d e k k e r , d e k n u t h ，m a m a l c o l m 幂n w k a h a n ，见参考文献 5 1 【2 2 】，【2 8 】， 【2 0 和 2 1 】。d m p r i e s t 和y k z h u 对于蒸馏算法( 与补偿算法类似) 做了细致分析，见【3 9 】， 【4 0 】， 5 3 1 着i 5 4 。美国加利福尼亚大学伯克利分校的j d e m m e l 教授也做了类似的研究， 见 7 仟口 8 】。最近德国汉堡大学r m r u m p 教授和日本早稻田大学t o g i t a 教授在前人基础 上对补偿算法做了系统的总结，并将称为无误差变换技术( e r r o r - f r e et r a n s f o r m a t i o n ) ，具 体见参见文献【4 3 】， 4 4 1 和【3 5 】。 在数值计算中，多项式函数值及其导数值的评估是一类经典问题，并且该问题的 大多数算法已经被深入研究。例如在误差分析领域著名的w i l k i n s o n 教授对评估多项式 的h o m e r 算法给出了向前和向后误差分析，见参考文献 5 0 讦口 5 1 】另外r h a m m e r 和n j h i g h a m 分别在文献【1 6 】的第四章和文献【1 8 】的第五章详细介绍了h o m e r 算法的误差分 析；r b a r r i o 在文献【3 】中对多项式评估给出了统一的误差界；i b m 的r e n e n k e l 提出了 特别的基于表格法( t a b l e b a s e dm e t h o d ) 评估多项式，见文献【1 0 】。近些年来法国居里大学 的e l a n g l o i s 教授和法国第六大学的s g r a i l l a t 昌1 教授根据r m r u m p 教授的理论提出了补 偿的h o r n e r 算法，算法精度能够达至i j i e e e 一7 5 4 双精度标准的最后一位有效位，与在双倍 双精度( d o u b l e - d o u b l e ) 标准下应用h o m e r 算法达到的精度相同但运算时问少、效率高，具 体参见文献【1 2 ，【1 3 】，【1 4 】，【2 3 】，【2 4 】及 2 5 1 。 多项式导数值的评估较函数值的评估复杂，常用算法如文献【1 8 】中所述。文献中 对算法的构成和舍入误差分析做了较为详细的论述。国际上对于这类问题的研究学者 有k h m f i l l e r 和o l v e r ，见参考文献 3 3 1 5 u 3 6 1 。但上述算法结果的相对误差小于相对条 件数和单位舍入误差量让( 见第二章) 的乘积，在条件数较大即问题病态的情况下很难达 第l 页 国防科学技术大学研究生院学位论文 到精度要求，这就需要一种更高精度的算法，我们在系统分析、借鉴前人科研成果的 基础上创造性的提出了一种新的补偿算法求解多项式函数的导数值，数值结果相对误 差小于相对条件数与o ( u 2 ) 乘积i j n _ k 单位舍入误差量t 。 h o m e r 算法通常用来计算幂指数基多项式函数值，对于c h e b y s h e v 基多项式相应的 有类似的算法一- - c l e n s h a w 算法，见参考文献【4 】。针对c l e n s h a w 算法，早期的研究学者 有j o l i v e r 和d e l l i o t t ，见参考文献【9 】和【3 7 】。最近，a s m o k t u n o w i c z 对c l e n s h a w 算法的向 后稳定性进行了分析，见【4 8 】。另外，t - h a s e g a w a 在m s k r z i p e k 提出的基于c l e n s h a w 算法 求多项式导数的算法基础上做了系统分析，见参考文献【1 7 】和【4 6 】。这些研究集中于分 析原算法的误差精度，而没有探索更高精度的算法，本文就是基于此要求并借鉴上述 学者的成果提出了补偿的c l e n s h a w 算法。 1 2 本文的工作及内容安排 本学位论文主要依据无误差变换技术提出高精度算法计算幂指数基多项式函数的 一阶导数值和c h e b y s h e v 基多项式的函数值。 第一章概述了国内外对于高精度计算理论发展情况，幂指数基多项式函数值和导 数值评估以及c h e b y s h e v 基多项式函数值计算等方面的研究现状，在此基础上提出了本 论文研究的方向和重点，同时给出了文章的内容安排。 第二章详细阐述了计算机浮点运算、舍入误差、i e e e 7 5 4 标准、条件数等相关概 念，介绍了s m r u m p 教授的误差自动转换技术发展情况，给出了基本加法和乘法的无 误差变换的算法，同时给出了至关重要的无误差变换恒等式定理，为后续算法精度分 析提供理论基础。 在第三章中，本文研究了一类基为幂指数的单变量多项式函数在一浮点值的函数 值及其导数值的精确评估文章首先系统分析了传统的评估多项式函数值及其导数值 的h o m e r 方法，然后基于无误差变换技术，提出了补偿的h o m e r 导数算法来提升传统的 求多项式一阶导数值的算法精度。该算法将传统算法计算步骤中的舍入误差记录下来， 通过分析得出传统算法最终数值结果的误差为三个分别以舍入误差及其组合为系数的 幂指数基多项式的和，进而将其补偿给传统算法数值结果，得到更精确的结果。最后 通过理论分析和数值实验证明了算法的高精度性和可靠性。 第四章在上一章的基础上研究了一类c h e b y s h e v 基的单变量多项式函数在一浮点值 的精确评估在简要介绍c h e b y s h e v 多项式和评估c h e b y s h e v 基多项式函数值的c l e n s h a w 算 法后，我们提出了类c h e b y s h e v 多项式和评估类c h e b y s h e v 基多项式函数值的类c l e n s h a w 算 第2 页 国防科学技术大学研究牛院学位论文 法的概念，为后续证明提供了理论基础。通过分析c l e n s h a w 算法计算步骤中舍入误差 的构成和累积，本文应用无误差变化技术将这些舍入误差提出，并以此为系数构建一 个c h e b y s h e v 基多项式，然后在理论上证明所构建的多项式的函数值即为原c l e n s h a w 算 法的数值误差。将新构建的c h e b y s h e v 基多项式的数值结果补偿回原数值结果，即为本 文提出的补偿c l e n s h a w 算法。误差理论分析和数值结果都展示了该算法的高精度性和 可靠性。 第3 页 国防科学技术大学研究生院学位论文 第二章预备知识 2 1 计算机浮点运算和舍入误差 绝大部分计算机内部以二进制方式工作而以十进制方式与用户进行交流，而转化 过程由计算机执行。但是由于计算机浮点数中尾数的位数是有限的，在进行科学计算 的时候，对一般实数总是要按舍入原则表示为浮点数，这就牵涉了一些小的舍入误差， 例如，o 1 这样的简单的实数也不能在二进制计算机中精确存储，通常我们不会注意这 一微小的转换误差，因为默认格式显示出来的还是0 1 。除了对输入数据进行舍入外， 在浮点数之间的四则运算结果同样也常常需要舍入。这些微小的舍入误差看似无关紧 要，但在大规模科学计算或病态问题中这些舍入误差可能积累到使计算结果失去意义 的地步，因此，有必要对浮点数算术运算的舍入误差进行系统的理论分析，并以此构 造更加精确可靠的数值算法。 首先给出浮点数系统的定义： 定义2 1 1f 表示浮点数全体的集合，它是实数的子集，即 f = o ) i - j ，：，= ：k m 矿q ) ( 2 1 1 ) 其中p 是计算机浮点系统的基数，当p = 2 ，8 ，1 0 ，1 6 时，分别表示二进制，八进制，十进 制和十六进制数；t 是字长表示精度；e 是指数，即阶数，有e m i n e e m a x ；尾数m 是 一个整数满足o mss t 一1 为了保证对任意浮点数，f 有唯一的表示形式，一般 采用规格化的科学记数法，这样一来就假设m 一1 。在f d g 非零的浮点数的范围 是伊m 饥一1 i f l 伊( 1 一p 一。) 这个界限就是浮点系统的下溢和上溢阀值。 我们根据定义2 1 1 可推出其改进形式 t f = o ) u ，：，= 士矿也p “ ( 2 1 2 ) 七= 1 其中整数d 七满足0 d 七sp l ，并且由于规格化，d 1 0 d x 为最高有效位，d t 为最低有 效位。 在给出了浮点系统的一般概念后，我们简介一下i e e e 7 5 4 标准一一这一国际上应用 最为广泛的浮点系统标准。i e e e 7 5 4 标准于1 9 8 5 年出版，定义了一个二进制浮点算术系 统，给出了如浮点数格式，基本浮点运算结果，舍入模式，规格化与非规格化数，保护位， 不同算术格式之间的转化等概念，详细资料见参考文献【1 9 】。形旺7 5 4 标准的浮点单精 第4 页 同防科学技术大学研究生院学位论文 度和双精度运算的参数在表2 1 中给出，其中：是有效位，s i z e 表示比特位( b i t s ) ，e p s i l o n 是 一个与单位相比可以忽略的整数( 使得g + 1 1 的最小浮点数e m ) ，h u g e 是最大的数( 上 溢阀值) ，t i n y 是最小的正数( 下溢阀值) ，p r e c i s i o n 是十进制精度，u = 1 2 5 1 一。表示单位 舍入误差( u n i tr o u n d o f f ) ，其具体说明见定理2 1 1 ，其它参数前文已经交代。 表2 1i e e e 7 5 4 标准单双精度浮点运算参数 m e e 7 5 4侈s i z e t e p s i l o np r e c i s i o n e m a xe m i n h u g et i n y u s i n g l e 2 3 2 2 42 2 3 61 2 81 2 52 1 2 82 1 2 65 9 6 1 0 8 d o u b l e 26 45 32 5 2 1 51 0 2 4 1 0 2 1 2 1 0 2 42 1 0 2 2 1 1 l 1 0 1 6 i e e e 7 5 4 标准的单精度和双精度格式的存储模式如图2 1 所示：其中8 为符号位，e 是 指数，m 是尾数，参见定义2 1 1 s i n g l e lb i t8 b i t s2 3 b i t s d o u b l e l b i tll b i t s5 2 b i t s 图2 1i e e e 7 5 4 标准单精度( 3 2 位) 和双精度( 6 4 位) 的存储形式 人们在进行高精度计算、处理极度病态问题或需要可靠计算结果的时候需要利用 t g r e e e 单精度或双精度更高的精度。故i e e e 定义了一个称为双倍扩充的更高精度，但 并不是所有的语言和编译器都能允许双倍扩充精度变量来计算和显示，所以工程中 一般采用应用软件来模拟更准确的运算，如下一节将要介绍的无误差变换( e r r o r - f r e e t r a n s f o r m a t i o n ) 。 下面介绍计算机的舍入误差。默认的舍入模式是舍入到最接近的数，即选择左右 较接近的那个浮点数。当距离相同时，采用舍入到偶数，即若位于左右两个浮点数中间 时，选择偶数的浮点数。应用这一模式的最大误差是最低有效位上的半个单位。其它 的舍入模式还有直接舍入，如向0 舍入( 截断) ，向+ 。o 舍入和向一。舍入。下面的定理 给出了在就近舍入模式下某一实数和其在浮点系统中的浮点数近似值的相对误差界。 定理的详细证明过程见参考文献【1 8 t h e o r e m2 2 。 第5 页 国防科学技术大学研究生院学位论文 定理2 1 1 若z i r 位于浮点数体统f 范围内，将实数z 的浮点数表示为似z ) ，则 ，f ( z ) = z ( 1 + 6 1 ) = 去，让，蚓u ( 2 1 3 ) 定理2 1 1 给出了实数z 和浮点数似z ) 的相对误差界限 i 华i - - - 州等胁 为了对数值算法进行舍入误差分析，科学计算领域通常应用以下基本浮点运算模 型： f l ( ao pb ) = ( ao p6 ) ( 1 + 6 ) ， 1 6 l 冬让，o p + ，一，x ， ( 2 1 4 ) 上述模型给出了因有限字长导致的计算机计算误差的过程，当ao p6 被计算和存储 时，计算机得到的是以舍入至l j f l ( ao p6 ) 机器字的形式来拟合ao pb ，然后将其存储， 并保证计算结果和真实值之间的相对误差界小于单位舍入误差t 本论文中还应用符 号o ，e ，o 表示浮点加、减和乘运算，如a0 6 = f l ( a + 6 ) ；同样应用符号z 加上标的形式 表示计算机浮点运算结果，如圣= 似z ) 。 在给出浮点数计算模型后，为了分析追踪模型中( 1 + 6 ) 因子对最终误差分析的影响， 我们给出了误差分析中常用的经典符号妒，具体见以下定理，详细说明见参考文献 1 8 1 l e m m a3 1 。 定理2 1 2 若民i t , p i = 4 - 1 ，对i = 1 ：且伽 1 ，则有 i i ( 1 + 魂) 店= 1 + ( 2 1 5 ) t = l 其中 1 0 1 尚= ： ( 2 1 6 ) 同样进一步方便分析误差传播，我们给出如下关系； ( 1 + 以) ( 1 + 岛) = 1 + p 七村， ( 2 1 7 ) 饥饥+ 1 ， ( 2 1 8 ) 饥 t m i n ( k ，j ) ， ( 2 1 9 ) 饥+ 竹+ 饥饥何 ( 2 1 1 0 ) 本论文中除了采用以上的一些符号表示方法还借鉴t s t e w a r t 4 7 提出的 符号： k = h ( 1 + 魂) 店，p i = 土1 ，1 6 以l 乱， ( 2 1 1 1 ) t = l 第6 页 国防科学技术大学研究生院学位论文 及相对关系式： = ， i k ： j 除了以上介绍的概念和符号外，在估计计算结果的误差过程中，我们还需一类至 关重要的概念，即扰动理论及条件数，在这一领域最为成功的莫过于j h w i l k i n s o n 教授 下面给出通用的条件数的定义，详见参考文献 4 1 】。 定义2 1 2 考虑求满足方程f ( z ，d ) = o 的z 值这类问题其中d 是一组解z 所依赖的数据， f 是关于z 和d 的函数，其中z 和d 可以是实数，向量，矩阵等。要求有问题有唯一解z ， 且关于依赖数据连续( w e l l p o s e d ) ，5 d 为输入数据扰动，妇为输出结果扰动，且有f ( z + 5 x ，d + 5 d ) = 0 。 相对条件毵 c o n d ，洲) - s u 舄黼黼( 2 1 1 2 5 d e o a a 1 ) dl | i i | j 绝对条件数： c o r t d a b s = 黜黼( 2 1 1 3 ) 其中d 是初始输入数据允许扰动的范围。 根据定义2 1 2 ，并结合本文研究对象的情况，我们给出简单的例子加以解释说明。 设，( z ) 是实变量z 的一个实值可微函数，要计算，( z ) 值，则即f ( d ，z ) = f ( z ) 一d = 0 ， 一个已知x 求d 的方程，这时根据定义2 1 2 ，并利用一个简单的对，的线性近似得到估 计f ( x + 5 x ) f ( x ) + 妇，7 ( z ) ，得至l j f 7 ( z ) i 为这一问题的绝对条件数，i ，7 ( z ) i i x l i f ( x ) l 为 相对条件数。 由以上条件数的分析，我们再引入优良算法具有的性质向后稳定性这一概念，见 参考文献【6 】，其定义如下： 定义2 1 3 设口冶( z ) 是含有舍入影响的，( z ) 的算法若对一切z 存在一个小的扰动6 z 使 得a l g ( x ) = f ( x + 如) ，则称。匆( z ) 是，( z ) 的向后稳定算法，妇称为向后误差。 由定义2 1 3 可以看出，若a l g ( ) 是向后稳定的，由于 e r r o r = i a l g ( x ) 一，( z ) i = i f ( z + 6 z ) 一，( z ) i i ，7 ( z ) i l ( i x l ( 2 1 1 4 ) 因为陋z h 总是很小，故当条件数i f 7 ( z ) i 不大时，误差将控制在很小的范围。我们提出的 算法应该尽可能都是向后稳定的算法。 第7 页 国防科学技术大学研究生院学位论文 2 2 无误差变换技术 无误差变换( e r r o rf r e et r a n s f o r m a t i o n ) 的思想很早就被提出，最近几年更是引起了学 术界广泛的注意。相当一部分计算浮点数向量和及内积的算法应用了这一技术。所谓 无误差变换简而言之就是： 定义2 2 1 设o + ，一，) ，口和6 是两个浮点数，且有z = f l ( a o6 ) ，这时存在 aob = z + y ( 2 2 1 ) 其中z 表示计算结果最好的浮点数近似，! ，表示精确的误差结果将浮点数对( n ，b ) 转化 为( z ，耖) ，且满足上式f 2 2 j j 即为无误差变换。 这一研究领域有早期的k n u t h 2 2 ，他提出了算法1 给出了两个浮点数的和的无误差变 换形式，算法需要6f l o p s ( 浮点运算量) 。d e k k e r 5 同样给出了一中带分支的浮点数和 的无误差转换算法，见算法2 ，算法需要4 f l o p s 。减法原理同加法是一样的。对于乘法首 先需要将输入参数分割成两部分，见算法3 及参考文献 5 】。其中r = r t l 2 1 ，为浮点系统 的尾数位。算法3 将一个浮点数a 分成两部分具有r 一1 个非零有效位的浮点数z 和y ，使 得a = z - i - y 。例如在i e e e 7 5 4 刃j 精度标准中，t = 5 3 ，r = 2 7 ，所以输出的数最多有2 6 位。 在此基础上v e l t k a m a p 【5 给出了算法4 。算法需要1 7 f l o p s 。 算法1 浮点数和的无误差变换形式 f u n c t i o n x , y = t w o s u m ( a x = f l ( a + b ) z = f l ( x - a ) y = f l ( a 一x 一乙) ) + ( b d 算法2 带分支的浮点数和的无误差变换形式 f u n c t i o n x , y = t w o s u m b r c ( a , 纠 x = f l ( a + b ) i f l a l l b i y = f l ( b - ( x 一口j j e l s e ，= l n 一( x - b ) ) 算法3 浮点数无误差变换形式分割算法 f u n c t i o n x , y = s p l i t ( a ) 第8 页 国防科学技术大学研究生院学位论文 f a c t e r = 2 + l z = f l ( f a c t o r 8 j x = f l ( z i z 口) ) y = f l ( a - x ) 算法4 浮点数乘积的无误差变换形式 f u n c t i o n x , y = t w o p r o d ( a ，纠 x = f l ( a ” 【口l ，a 2 】= s p l i t ( a ) 【6 l ，b 2 1 = s p l i t ( b ) y = f l ( a 2 b 2 一( ( ( z 一口l 6 1 ) 一0 2 6 1 ) 一口l 6 2 ) ) 下面定理给出了上述算法的性质，其在后续的证明中有重要意义。 定理2 2 1 设口，b f _ e x ，y f ，那么有 f j j b ，y 】= t w o s u m ( a ，6 ) ，z = f l ( a + 6 ) ，z + y = a + b ，y u l x l ，ysu l a + b | 例k ，引= t w o p r o d ( a ，6 ) ，z = f l ( a 6 ) ，z + y = a b ，y u k i ，y u l a b i 至此，关于浮点数的无误差变换的基本运算构建完成。在这一基本思想提出以 来，学者们掀起了长时间的研究高潮并取得了长足的发展，涌现了大量的研究论文和专 著，女1 ：i 3 9 1 ， 5 3 1 和【l 】，这些文章主要是集中讨论向量运算的无误差变换即所谓的蒸馏算 法( d i s t i l l a t i o n a l g o r i t h m ) 。近年来r u m p ，o g i t a 和o i s h i 教授系统总结了无误差变换的思想 并提出了向量和和向量内积的相应无误差变换形式，提出了达到任意精度的向量和算 法，极大的推动了无误差变换的理论发展和应用，主要内容参见 4 3 】， 4 4 】，【4 5 】和【3 5 】。 第9 页 国防科学技术大学研究生院学位论文 第三章幂指数基多项式函数值及其导函数值的高精度计算 本章主要研究了基于幂指数基多项式的函数值及其导数值的无误差变换技术，依 据h o m e r 算法构造高精度的评估算法，并从理论上分析了该类算法在问题不是过于病态 的情况下精度达到了最后一位有效位，数值实验同样证实了我们的结论。 3 1 多项式及其导数值的评估 这一节简述计算多项式函数值和导数值的h o m e r 方法，为后续章节提供理论基础， 详细介绍见参考文献 1 8 1 。 3 1 1 多项式函数值计算 经典的评估多项式 p ( x ) = a o + a l x + + 矿( 3 1 1 ) 函数值的算法是h o m e r 算法，如下是其详细规则： 算法5h o r n e r 算- 法 f u n c t i o nq o = h o r e n e r ( p ，z ) 锄( z ) = a n 如ri = n 一1 ：- 1 ：0 吼( z ) = x q i d - l ( x ) + o t i e n d p ( x ) = q o = q o ( x ) 为了分析h o m e r 算法的舍入误差，我们引入符号 ，见预备知识公式( 2 1 1 1 ) ，根据 上述算法得到 蠡= ( x q i - i - 1 ( 1 + & ) + ) ( 1 + ) = ( o 磊+ 1 ) + a i 其中i 文i ，例u ， = ( 1 + 文) ，啦= a d z ) 根据缸= = ，展开迭代公式得到 哦= o t n g g n - i ( 2 ( n i ) ) q - r 、佗t = 一i 1q z t 一( 2 ( 一t ) + 1 ) ( 3 1 2 ) 第1 0 页 国防科学技术大学研究生院学位论文 当i = o 时得到最后的计算值表达式： n 一1t ln 一1 0 0 = 【( 1 + 彭) ( 1 + 妨) ( 1 + 嘭) 】q t + 【n ( 1 + 如) ( 1 + 嘭) 】z n i = o j = o j = o ( 3 1 3 ) = ( 1 + 0 1 ) o 0 + ( 1 + 0 3 ) a l x + + ( 1 + 6 1 2 竹一1 ) 一l z 佗一1 + ( 1 + 0 2 n ) 扩 = o 0 + a l x + + o n 一1 x n 一1 + z n 根据定理2 1 2 及公式3 1 3 我们给出h o m e r 算法的误差界： n i p ( x ) 一0 0 l 协l a l l x l = 7 2 n p ( i x l ) ( 3 1 4 ) i = o 其中钿为应用h o m e r 算法得到的带有误差的数值结果，定义了新的多项式函数： 佗 痧( z ) = 蚓 ( 3 1 5 ) 同样算法的相对误差如下 警素掣协黜= 协啪坳i p ( z ) 一。l p ( z ) i 。“、，“7 其中 删舭，= 黜= 超嬲 为多项式求值的相对条件数。 3 1 2 多项式导数值计算 ( 3 1 6 ) ( 3 1 7 ) 根据上一节的多项式函数值的计算方法，我们显然可以通过将多项式( 3 1 1 ) 求导得 出多项式导数公式，然后根据表达式应用h o m e r 算法计算多项式导数值。但是h i g h i a m 在 参考文献 1 8 1 中给出了一个更加有效的方法： 首先定义一个新的多项式 q ( x ) = fq i x “ ( 3 1 8 ) t = 1 其中吼= 吼( 口) ，h o m e r 算法( 算法5 ) 的每一计算步骤生成。这时我们很容易得到 p ( x ) = ( z a ) q ( x ) + q o( 3 1 9 ) 这类似于合成多项式除法，显然我们有 p t ( 口) = q ( a ) 对口( z ) 应用同样的处理方法，我们的将会得到多项式p ( 口) 在值a 处的克阶导数值。下面给 出求多项式导数值的算法： 第l 1 页 国防科学技术大学研究牛院学位论文 算法6 求后阶导数的h o m e 得数算法 f u n c t i o nt = h o r e n e r d i f k ( p ，a ) t o = o t n ，孔= t 2 = = 死= 0 加，j = 佗一1 ：- 1 ：0 归ri = m i n ( k ，n j ) ：- 1 ：1 死= 8 正+ 正一1 e n d t o = a t o + o t i e n d m = 1 力r j = 2 ：k m = m 木j t j = m 木t j e n d 其中互是多项式第i 阶导数值。对于一阶导数的特殊情况我们简化算法为： 算法7 求一阶导数的h o m e 厂导数算法 f u n c t i o ny o = h o r e n e r d i f ( p ，a ) q n = c z n y n = 0 如，i = 铊一1 ：- 1 ：0 玑= a y + l + q i + l 口t = a q i + l + o l i e n d 其中o $ l q o 分别为多项式在a 点处的导数值和函数值。显然上述算法可以分解成两部分： 一部分是对多项式p ( z ) 应用h o m e r 算法得到一系列数吼( z ) ，记为多项式g ( z ) 的系数，然后 对这一多项式再次应用h o m e r 算法得到多项式p ( z ) 在点a 处的一阶导数值。 对于一阶导数的h o m e r 算法，有 i ( 口) 一如i 2 佗乱( + o ( u 2 )( 3 1 1 0 ) 其中舶是应用算法7 得到多项式一阶导数的带有误差数值结果，f i ( a ) = 廷1 i q 知叫七 详细的推导过程见参考文献【1 8 】同样对于多项式的后阶导数可以推出类似的误差界。 第1 2 页 同防科学技术大学研究生院学位论文 3 2 评估多项式值的补偿算法 对于评估多项式函数值的h o m e r 算法，公式( 3 1 6 ) 给出了计算结果的相对误差界， 可见计算结果的精度小于机器提供的精度。一般情况下条件数c o n d ( p ，z ) 不大时，对结 果精度要求不高时影响不大；但在某些对精度要求较高或问题病态情况下，如多项 式函数p ( z ) 在z 处病态( 在多项式重根的附近) ，误差界 2 n c o n d ( , p ，z ) 可能会大于1 ，计 算结果失去意义，这时高精度的算法就很有必要。应用较为广泛是双倍精度( d o u b l e d o u b l e ) 算法，见参考文献1 2 6 ，其核心思想是模拟两倍的i e e e 双精度，算法可以应用于计 算方法的多个领域。最近s g r a i l l a t 禾i p l a n g l o i s 针对多项式评估提出的一种补偿h o m e r 算 法( c o m p e n s a t e dh o m e ra l g o r i t h m ) ，参见【1 2 】算法应用无误差变换技术，简单可靠，精度 高，速度快，不失为一种新的计算多项式函数值的可靠算法。 与应用双倍精度( d o u b l e d o u b l ea r i t h m e t i c ) 的h o m e r 算法达到的精度一样，补偿h o r n e r 算 法( c o m p h o m e r ) 的相对误差满足： i c o m p h o r i n e 两r ( p r , x ) - 一p ( x ) l 冬u + c o n d ( p , x ) o ( u 2 ) ( 3 2