曹 靖 ，李建平

（1.中國(guó)科學(xué)院a.大氣物理研究所，b.研究生院，北京100029；2.天津理工大學(xué)理學(xué)院，天津300384）

“大數(shù)吃小數(shù)”現(xiàn)象[1]是數(shù)值運(yùn)算中常見(jiàn)的一種影響計(jì)算精度的現(xiàn)象，當(dāng)以計(jì)算機(jī)計(jì)算一個(gè)實(shí)數(shù)a與實(shí)數(shù)b≠0的代數(shù)和時(shí)，如果|b|相對(duì)于|a|小到一定程度，會(huì)出現(xiàn)a+b=a的現(xiàn)象，一般稱(chēng)作數(shù)b被數(shù)a“吃掉”了.“大數(shù)吃小數(shù)”現(xiàn)象一個(gè)典型的例子[2]就是計(jì)算N個(gè)實(shí)數(shù)的累加和3，…，N.若以正常順序累加，則自a2之后所有元素均被a“1吃掉”，累加結(jié)果為如果N取值很大，那么很可能是個(gè)較大的數(shù)，這樣數(shù)值計(jì)算就會(huì)有一個(gè)較大的誤差.另一個(gè)例子是在數(shù)值運(yùn)算的網(wǎng)格剖分中，如果剖分選取的分辨率（步長(zhǎng)）過(guò)小，那么計(jì)算過(guò)程中它有可能會(huì)被“吃掉”，出現(xiàn)網(wǎng)格的某些節(jié)點(diǎn)無(wú)法生成的情況，從而進(jìn)一步影響數(shù)值運(yùn)算精度.

目前，已存在一些避免“大數(shù)吃小數(shù)”現(xiàn)象的方法[2，3].做N個(gè)實(shí)數(shù)的累加和時(shí)，可調(diào)換相加順序，將絕對(duì)值較小的數(shù)先進(jìn)行累加，以避免“大數(shù)吃小數(shù)”現(xiàn)象，但仍存在2種風(fēng)險(xiǎn)：1.前面“小數(shù)”的累加值相比于后面“大數(shù)”仍然足夠小到被“吃掉”，調(diào)換順序起不到作用；2.前面若干“小數(shù)”的累加值遠(yuǎn)大于后面的“大數(shù)”，從而將其“吃掉”，使運(yùn)算結(jié)果更加不可預(yù)料.可見(jiàn)，僅靠目前的方法并不能完全避免危及精度的現(xiàn)象發(fā)生.為真正避免此類(lèi)現(xiàn)象的發(fā)生，找到準(zhǔn)確判斷其發(fā)生的界限，就顯得十分迫切和必要.本研究針對(duì)兩同符號(hào)數(shù)相加的問(wèn)題，通過(guò)數(shù)值試驗(yàn)與理論分析，給出“大數(shù)吃小數(shù)”現(xiàn)象發(fā)生的界限，從而為實(shí)際運(yùn)算提供理論指導(dǎo).

1 機(jī)器雙精度下“大數(shù)吃小數(shù)”的界限

1.1 兩正數(shù)相加的數(shù)值試驗(yàn)

設(shè)A＞0為一“大數(shù)”，CA＞0為可被A“吃掉”的最大“小數(shù)”.在計(jì)算機(jī)雙精度下，分別應(yīng)用Matlab、Fortran以及C等3種計(jì)算機(jī)語(yǔ)言進(jìn)行大量數(shù)值試驗(yàn)，發(fā)現(xiàn)3種語(yǔ)言計(jì)算出的A與CA的關(guān)系都呈現(xiàn)出相同的規(guī)律.以Matlab為例，給出A∈[0.1，10]時(shí)的部分試驗(yàn)結(jié)果，如圖1所示.由圖1可見(jiàn)，A所處的區(qū)域可被分為若干區(qū)間段，在每個(gè)區(qū)間段內(nèi)，雖然A的取值不同，但被其吃掉的最大“小數(shù)”CA卻十分接近，從而圖形呈現(xiàn)階梯形狀.特別地，在每個(gè)區(qū)間段內(nèi)，lg A與lg（CA/A）之間呈現(xiàn)線性關(guān)系，因此，下面考慮應(yīng)用這種分段線性關(guān)系，推導(dǎo)出CA的計(jì)算公式.

圖1 機(jī)器雙精度下A與CA及l(fā)g A與lg（CA/A）的關(guān)系Fig.1 A versus CAand lg A versus lg（CA/A），under double machine precision

1.2 兩正數(shù)相加lg（CA/A）關(guān)于lg A的分段線性擬合

為考察lg（CA/A）關(guān)于lg A的分段線性關(guān)系，進(jìn)行了大量數(shù)值試驗(yàn).首先，根據(jù)試驗(yàn)數(shù)據(jù)計(jì)算出lg（CA/A）關(guān)于lg A斜率的變化規(guī)律，將lg A所屬區(qū)間進(jìn)行分段，即估計(jì)出試驗(yàn)所涉及的每段線性區(qū)間的左、右頂點(diǎn)坐標(biāo).然后，分別計(jì)算出lg A每個(gè)分段區(qū)間的長(zhǎng)度，以及區(qū)間內(nèi)lg（CA/A）關(guān)于lg A線性擬合的斜率.部分試驗(yàn)結(jié)果見(jiàn)表1.

表 1 機(jī)器雙精度下，當(dāng) A∈（10-2，102）時(shí)，lg（CA/A）與lg A線性關(guān)系Tab.1 Linear relationship between lg（CA/A）and lg A，when A∈（10-2，102），under double machine precision

由表1可見(jiàn)，不同區(qū)間段之間滿足十分相似的線性關(guān)系：lg A分段區(qū)間長(zhǎng)度都近似為0.301，并且，每個(gè)區(qū)間段內(nèi)，lg（CA/A）關(guān)于lg A線性擬合的斜率都接近于-1，且 lg（CA/A）均在[-16.256，-15.955]范圍內(nèi)單調(diào)遞減.由此便可估計(jì)lg（CA/A）關(guān)于lg A的分段線性擬合函數(shù).

首先，給出每個(gè)分段線性區(qū)間左頂點(diǎn)坐標(biāo).由表1可見(jiàn)，lg A=0為某一分段區(qū)間的左頂點(diǎn)，又因每個(gè)區(qū)間段長(zhǎng)度均近似為0.301，則可記“大數(shù)”A屬于第k（k∈Z）段區(qū)間，其中 k=floor（lg A/0.301），floor表示向負(fù)無(wú)窮方向取整.并且，此第k段區(qū)間的左頂點(diǎn)坐標(biāo)應(yīng)為（0.301k，-15.955），k∈Z.然后，以表示 CA的分段線性擬合近似值，若A屬于第k段區(qū)間，結(jié)合上述試驗(yàn)中區(qū)間段內(nèi)lg（CA/A）關(guān)于lg A線性擬合的斜率以及左頂點(diǎn)坐標(biāo)，得此區(qū)間段內(nèi)lg（/A）關(guān)于lg A的線性擬合函數(shù)為

為驗(yàn)證式（1）的擬合效果，隨機(jī)選取不同的“大數(shù)”A，將數(shù)值試驗(yàn)得到的CA值與由式（1）計(jì)算出的擬合值進(jìn)行比較，見(jiàn)表2.表2數(shù)據(jù)說(shuō)明式（1）具有良好的擬合效果.

表2 機(jī)器雙精度下與CA比較Tab.2 Comparison between and CAfor random，under double machine precision

表2 機(jī)器雙精度下與CA比較Tab.2 Comparison between and CAfor random，under double machine precision

A數(shù)值試驗(yàn)結(jié)果CA 擬合值C A 1.00 1.109×10-16 1.109×10-16 16.13 1.776×10-15 1.774×10-15 282.60 2.839×10-14 2.838×10-14 3 596.00 2.269×10-13 2.270×10-13 7.92×106 4.653×10-10 4.645×10-10 1.00×10-2 8.670×10-19 8.670×10-19 3.45×10-3 2.167×10-19 2.168×10-19 5.88×10-5 3.383×10-21 3.388×10-21

1.3 兩負(fù)數(shù)相加情況

設(shè)A′＜0為一“大數(shù)”，應(yīng)尋找C′A為全體負(fù)數(shù)中可被A′吃掉的最小的“小數(shù)”，同時(shí)也是絕對(duì)值最大的“小數(shù)”.大量數(shù)值試驗(yàn)結(jié)果表明，當(dāng)|A′|=A時(shí)，總有|C′A|=CA.因此，可將上述兩正數(shù)相加情況的分段線性擬合結(jié)果直接推廣至兩負(fù)數(shù)相加情況.令′A為C′A的擬合值，則結(jié)合式（1）可得

1.4 雙精度下“大數(shù)吃小數(shù)”界限的近似公式

結(jié)合式（1）與式（2），可得機(jī)器雙精度下兩同號(hào)數(shù)相加時(shí)，“大數(shù)吃小數(shù)”現(xiàn)象發(fā)生界限的統(tǒng)一近似公式.

結(jié)論1 在機(jī)器雙精度下進(jìn)行運(yùn)算，對(duì)于任意符號(hào)相同的兩實(shí)數(shù)A與B，當(dāng)且僅當(dāng)|B|≤CA時(shí)，機(jī)器運(yùn)算結(jié)果為A+B=A，其中

2 機(jī)器單精度下“大數(shù)吃小數(shù)”的界限

經(jīng)與雙精度類(lèi)似的過(guò)程分析，可得機(jī)器單精度下“大數(shù)吃小數(shù)”界限的統(tǒng)一近似公式.

結(jié)論2 在機(jī)器單精度下進(jìn)行運(yùn)算，對(duì)于任意符號(hào)相同的兩實(shí)數(shù)a與b，當(dāng)且僅當(dāng)|b|≤ca時(shí)，機(jī)器運(yùn)算結(jié)果為a+b=a，其中

表3給出對(duì)于隨機(jī)選出的不同“大數(shù)”a，可被其“吃掉”的最大“小數(shù)”ca的數(shù)值試驗(yàn)值與由式（4）計(jì)算出的近似值ca的比較結(jié)果，這些結(jié)果說(shuō)明式（4）有良好的估計(jì)效果.

表3 機(jī)器單精度下ca試驗(yàn)值與式（4）近似值ca比較結(jié)果Tab.3 Comparison between caand cafor random，under single machine precision

3 任意機(jī)器精度“大數(shù)吃小數(shù)”界限的普適近似公式

公式（3）與公式（4）非常相似，只有指數(shù)上的一個(gè)常系數(shù)不同，單精度為7.225，雙精度為15.955.若記所選機(jī)器精度的二進(jìn)制有效位數(shù)為n，則機(jī)器單雙精度下n分別為24和53[4].注意到24lg 2≈7.224 7，53lg 2≈15.954 6，也就是說(shuō)，此常系數(shù)可由n近似推算出來(lái)，為nlg 2，由此推測(cè)CA與ca的近似表達(dá)式與n有關(guān).因此，在具有n位二進(jìn)制有效數(shù)字的機(jī)器下進(jìn)行運(yùn)算，將被“大數(shù)”a“吃掉”的“小數(shù)”界限記為ca，n，結(jié)合結(jié)論1和結(jié)論2，得任意機(jī)器精度下“大數(shù)吃小數(shù)”界限的普適近似公式.

結(jié)論3 在具有n位二進(jìn)制有效數(shù)字的機(jī)器精度下進(jìn)行運(yùn)算，對(duì)于任意符號(hào)相同的兩實(shí)數(shù)a與b，當(dāng)且僅當(dāng)|b|≤ca，n時(shí)，機(jī)器運(yùn)算結(jié)果為 a+b=a，其中

4 結(jié)論

通過(guò)數(shù)值試驗(yàn)，為數(shù)值運(yùn)算中避免“大數(shù)吃小數(shù)”現(xiàn)象給出了界限，結(jié)論3適用于任意機(jī)器精度和任意計(jì)算機(jī)語(yǔ)言，可方便應(yīng)用于實(shí)際的數(shù)值計(jì)算中.下一步考慮利用二進(jìn)制對(duì)位相加過(guò)程[5]，對(duì)于“大數(shù)吃小數(shù)”問(wèn)題進(jìn)行理論分析，并給出更為嚴(yán)格的理論界限.

另外，關(guān)于“大數(shù)吃小數(shù)”現(xiàn)象，本文只討論了同號(hào)的情況，可應(yīng)用類(lèi)似的方法研究?jī)僧愄?hào)數(shù)相加的情況，找到“大數(shù)吃小數(shù)”現(xiàn)象發(fā)生的統(tǒng)一界限公式.

[1]劉目樓，汪卉琴.數(shù)值分析[M].北京：冶金工業(yè)出版社，2005.

[2]李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].北京：清華大學(xué)出版社，2008.

[3]徐士良.計(jì)算機(jī)常用算法[M].北京：清華大學(xué)出版社，2005.

[4]朱亞超.基于IEEE754的浮點(diǎn)數(shù)存儲(chǔ)格式分析研究[J].計(jì)算機(jī)與信息技術(shù)，2006（9）：50-52.

[5]GOLDBERGD.ComputerArithmetic[M].Amsterdam：ElsevierScience，2003.