同符號數(shù)相加“大數(shù)吃小數(shù)”的界限：數(shù)值試驗

曹 靖 ，李建平

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

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

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

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

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

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

圖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ù)值試驗.首先，根據(jù)試驗數(shù)據(jù)計算出lg（CA/A）關(guān)于lg A斜率的變化規(guī)律，將lg A所屬區(qū)間進(jìn)行分段，即估計出試驗所涉及的每段線性區(qū)間的左、右頂點坐標(biāo).然后，分別計算出lg A每個分段區(qū)間的長度，以及區(qū)間內(nèi)lg（CA/A）關(guān)于lg A線性擬合的斜率.部分試驗結(jié)果見表1.

表 1 機(jī)器雙精度下，當(dāng) A∈（10-2，102）時，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可見，不同區(qū)間段之間滿足十分相似的線性關(guān)系：lg A分段區(qū)間長度都近似為0.301，并且，每個區(qū)間段內(nèi)，lg（CA/A）關(guān)于lg A線性擬合的斜率都接近于-1，且 lg（CA/A）均在[-16.256，-15.955]范圍內(nèi)單調(diào)遞減.由此便可估計lg（CA/A）關(guān)于lg A的分段線性擬合函數(shù).

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

為驗證式（1）的擬合效果，隨機(jī)選取不同的“大數(shù)”A，將數(shù)值試驗得到的CA值與由式（1）計算出的擬合值進(jìn)行比較，見表2.表2數(shù)據(jù)說明式（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ù)值試驗結(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ù)”.大量數(shù)值試驗結(jié)果表明，當(dāng)|A′|=A時，總有|C′A|=CA.因此，可將上述兩正數(shù)相加情況的分段線性擬合結(jié)果直接推廣至兩負(fù)數(shù)相加情況.令′A為C′A的擬合值，則結(jié)合式（1）可得

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

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

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

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

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

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

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

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

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

公式（3）與公式（4）非常相似，只有指數(shù)上的一個常系數(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，也就是說，此常系數(shù)可由n近似推算出來，為nlg 2，由此推測CA與ca的近似表達(dá)式與n有關(guān).因此，在具有n位二進(jìn)制有效數(shù)字的機(jī)器下進(jìn)行運算，將被“大數(shù)”a“吃掉”的“小數(shù)”界限記為ca，n，結(jié)合結(jié)論1和結(jié)論2，得任意機(jī)器精度下“大數(shù)吃小數(shù)”界限的普適近似公式.

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

4 結(jié)論

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

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

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

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

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

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

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