基于SCILAB的多精度算法研究與實現(xiàn) *

蘭 靜，劉文超，姜 浩,林文強

(1.重慶工商大學(xué)融智學(xué)院，重慶 404100；2.國防科技大學(xué)計算機(jī)學(xué)院，湖南 長沙 410073)

1 引言

高性能處理器被廣泛應(yīng)用于信息化聯(lián)合作戰(zhàn)、航空航天、模擬核試驗、氣象預(yù)報和搶險救災(zāi)等領(lǐng)域。近年來，我國自主研發(fā)的處理器發(fā)展迅速，這對于我國的國家安全具有重要意義，有利于完善我國經(jīng)濟(jì)信息系統(tǒng)以及國防領(lǐng)域的安全和保密問題。隨著計算機(jī)技術(shù)的不斷發(fā)展，高端計算已經(jīng)開始從單一追求高性能轉(zhuǎn)變?yōu)橹铝τ诰C合高效能的提升，以解決當(dāng)前高性能計算領(lǐng)域所面臨的實用性能和可編程性問題。浮點數(shù)運算對于計算機(jī)性能影響顯著，在浮點運算過程中存在舍入誤差，計算規(guī)模較小時，這些誤差不會造成明顯影響?，F(xiàn)代科學(xué)工程計算一般涉及到大規(guī)模計算，在這種大規(guī)模、長時程的數(shù)值計算中，由于舍入誤差的累積效應(yīng)存在，可能導(dǎo)致數(shù)值計算結(jié)果不準(zhǔn)確、不可靠，甚至得到完全錯誤的結(jié)果。高精度算法庫的設(shè)計和實現(xiàn)，是解決大規(guī)模數(shù)值計算的可靠性和穩(wěn)定性，以及提升國產(chǎn)處理器性能的有效途徑之一。

設(shè)計和實現(xiàn)高精度的處理器，可以從硬件和軟件2個方面進(jìn)行。在20世紀(jì)60年代到80年代，IBM公司研發(fā)出一系列處理器，最大特點是開發(fā)出了十六進(jìn)制的四精度浮點運算硬件邏輯，但是采用的并不是IEEE-754 浮點算術(shù)標(biāo)準(zhǔn)，而是IBM獨立的標(biāo)準(zhǔn)，其兼容性和推廣性有待提高。到了1998年，IBM公司研制的S/390G5處理器[1]正式采用了IEEE-754標(biāo)準(zhǔn)。隨后，IBM的z結(jié)構(gòu)處理器得到了硬件支持[2]。軟件實現(xiàn)高精度運算主要分為2大類：多部分(Multiple-Component)模式和多數(shù)字(Multiple-Digit)模式。多部分模式把一個高精度數(shù)據(jù)分為若干個標(biāo)準(zhǔn)的浮點數(shù)，其數(shù)值由這若干個浮點數(shù)的數(shù)值求和得到。每個浮點數(shù)擁有自己獨立的指數(shù)和尾數(shù)。多部分模式數(shù)據(jù)的精度是工作精度的整數(shù)倍，數(shù)據(jù)之間的運算可以通過浮點操作模擬實現(xiàn)。因此，多部分模式高精度運算具有良好的可移植性，而且可以充分利用高性能處理器所提供的浮點性能，計算速度相對于多數(shù)字模式高精度運算有明顯提高。本文主要應(yīng)用多部分模式設(shè)計和實現(xiàn)高精度算法庫。

自主可控是國家信息化建設(shè)的關(guān)鍵環(huán)節(jié)，是保護(hù)信息安全、國防安全的重要目標(biāo)之一。面向自主可控處理器平臺的軟件開發(fā)是急需解決的問題。Matlab是應(yīng)用較為廣泛的數(shù)值分析和圖形圖像處理商用軟件，為廣大的科研工作者所熟知和使用。一方面，在某些涉及國家安全的核心科研領(lǐng)域，如核武器模擬、大型航天器設(shè)計等，自主可控軟件是首要要求；另一方面，Matlab的內(nèi)部核心代碼不公開，在國產(chǎn)自主處理器上的移植也需要授權(quán)并由該公司主導(dǎo)開發(fā)?？紤]到SCILAB是非常接近于Matlab的開源軟件，在面向自主可控國產(chǎn)處理器平臺上，基于SCILAB進(jìn)行開發(fā)是首要選擇。SCILAB是法國國立信息與自動化研究院的科學(xué)家開發(fā)的開源科學(xué)計算軟件。此軟件的使用、修改和分發(fā)不受許可證的限制，自主可控。這樣就可以依靠自身研發(fā)設(shè)計，全面掌握產(chǎn)品核心技術(shù)，實現(xiàn)信息系統(tǒng)從硬件到軟件的自主研發(fā)、生產(chǎn)、升級、維護(hù)的全程可控。

基于此，本文利用無誤差變換技術(shù)，采用double-double數(shù)據(jù)格式，以SCILAB運算符重載和函數(shù)重載為手段，設(shè)計和實現(xiàn)面向SCILAB開源軟件的多精度的數(shù)值算法庫。應(yīng)用本文設(shè)計的多精度算法庫，科研工作者可以進(jìn)行更高精度的數(shù)值計算和模擬，對數(shù)值計算的舍入誤差進(jìn)行有效控制，從而使得計算結(jié)果的可靠性和穩(wěn)定性得到進(jìn)一步的保障。該算法庫主要面向國產(chǎn)自主設(shè)計的處理器，具有較好的可移植性，可以方便地將Matlab的數(shù)值模擬程序移植到國產(chǎn)處理器平臺。本文是文獻(xiàn)[3]的擴(kuò)展版本。

2 多精度算法研究與實現(xiàn)

2.1 無誤差變換

2005年，日本學(xué)者Ogita和Oishi以及德國學(xué)者Rump[4]在前人工作的基礎(chǔ)上提出了無誤差變換這一概念。設(shè)°∈{+,-,×}，2個浮點數(shù)(a,b)∈F且有x=fl(a°b)∈F，其中F為浮點數(shù)集合。在沒有上下溢出,舍入模式為就近舍入時，存在：

(a°b)=x+y

(1)

其中,x表示浮點計算結(jié)果，y∈F表示精確的舍入誤差結(jié)果。將浮點數(shù)對(a,b)轉(zhuǎn)化為浮點數(shù)對(x,y)，且滿足式(1)的過程即為浮點數(shù)加、減和乘法運算的無誤差變換。

(1)2個浮點數(shù)加法的無誤差變換[5]如算法1所示。

算法1FastTwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a⊕b；

步驟2y=b?(x?a)。

(2)兩個浮點數(shù)乘法的無誤差變換[6]如算法2所示。

算法2TwoProd

輸入：a,b。

輸出：x,y。

步驟1x=a?b；

步驟2[a1,a2]=Split(a)；

步驟3[b1,b2]=Split(b)；

步驟4y=a2?b2-(((x-a1?b1)-a2?b1)-a1?b2)。

其中，Split表示把1個浮點數(shù)拆成1對浮點數(shù)。

2.2 double-double算術(shù)

由美國勞倫斯伯克利國家重點實驗室的Hida等人[7]開發(fā)的軟件庫 QD(Quotient-Difference algorithm)應(yīng)用非常廣泛，該軟件首次提出了雙倍雙精度(double-double)和四倍雙精度(quad-double)格式。日本理化研究院的Nakata[8]開發(fā)的軟件Mpack針對LAPACK(Linear Algebra PACKag)和BLAS(Basica Liner Algebra Subprograms)實現(xiàn)了高精度的MLAPACK(Multiple precision arithmetic of LAPACK)和MBLAS(Multiple Precision arithmetic of BLAS)，其部分代碼是基于QD高精度軟件庫設(shè)計的。2010年，香港科技大學(xué)的Lu等人[9]在GPU硬件加速平臺上基于QD設(shè)計了GQD(GPU QD)高精度函數(shù)庫。同年，日本筑波大學(xué)的Mukunoki等人[10]在GPU上實現(xiàn)了BLAS函數(shù)庫四倍精度版本。國防科技大學(xué)的姜浩等研究人員[11-14]設(shè)計實現(xiàn)了多項式求根的高精度補償數(shù)值迭代法。double-double算法設(shè)計簡單，易于改進(jìn)當(dāng)前的數(shù)值算法，算法運行較快，適用于本課題的研究。

假定xh和xl是2個浮點數(shù)(xh,xl∈F)，x表示xh和xl的沒有舍入操作的精確和，即:

x=xh+xl

其中，xh和xl沒有交疊，且xh和xl是符合IEEE-754標(biāo)準(zhǔn)的雙精度double格式數(shù)，那么x就是double-double格式數(shù)。

2.3 運算符重載

運算符重載，就是對已有的運算符進(jìn)行重新定義，賦予其另外一種功能，以滿足不同的數(shù)據(jù)類型。重新定義的運算符重載函數(shù)的作用與內(nèi)置運算符的作用類似，下面重點介紹SCILAB中的重載實現(xiàn)規(guī)則。

二元操作數(shù)變換規(guī)則：

%〈first_operand_type〉_〈op_code〉_〈second_operand_type〉

一元操作數(shù)變換規(guī)則：

%〈operand_type〉_〈op_code〉

其中，operand_type表示操作數(shù)的數(shù)據(jù)類型,op_code表示操作符(如+，-)的類型。表1和表2列出了部分對應(yīng)規(guī)則。

本文將對表3中所列出的運算符進(jìn)行重載。

根據(jù)運算符重載的規(guī)則，以加法重載的實現(xiàn)為

Table 1 Cross reference list to operand_type

Table 2 Cross reference list to op_code

Table 3 Overloading operational character

例，對運算符重載的過程進(jìn)行展示。把浮點數(shù)a，b分割成高位H和低位L，分別記為ah，al，bh，bl，然后實現(xiàn)double-double的加法。在此基礎(chǔ)上使用重載命令對加法進(jìn)行重載?；赿ouble-double數(shù)據(jù)格式的加法重載過程為：

functionx=add_dd_dd(ah,al,bh,bl)

[s1,s2]=TwoSum(ah,bh);

[t1,t2]=TwoSum(al,bl);

s2=s2+t1;

t1=s1+s2;

s2=s2-(t1-s1);

t2=t2+s2;

[rh,rl]=FastTwoSum(t1,t2);

x=list(rh,rl);

endfunction

deff(′x=%l_a_l(a,b)′,′x=add_dd_dd(a(1),a(2),b(1),b(2))′)

其中，deff表示重載操作。

2.4 函數(shù)重載與構(gòu)建

SCILAB中已經(jīng)定義了如求絕對值(abs)、求平方根(sqrt)等函數(shù)，但是這些函數(shù)只適用于普通的浮點數(shù)運算。本文將重新定義這些函數(shù)，使其適用于雙倍浮點數(shù)精度的運算?；赿ouble-double數(shù)據(jù)格式，對表4中的函數(shù)實現(xiàn)了重載。

Table 4 Overloading functions

以equal函數(shù)為例，本文基于double-double格式對其進(jìn)行重載，重載具體實現(xiàn)過程為：

functiony=equal_dd(a,b)

y=(a(1)==b(1))&(a(2)==b(2));

endfunction

deff(′y=equal(a,b)′,′y=equal_dd(a,b′));

其中,equal_dd是基于double-double格式重新定義的函數(shù)，使用deff進(jìn)行重載操作。

3 數(shù)值實驗與分析

本節(jié)中，所有的數(shù)值實驗都是在 IEEE-754 標(biāo)準(zhǔn)[15]雙精度下進(jìn)行的，用在數(shù)值實驗中的多項式的系數(shù)和估值點都是浮點數(shù)的形式。精度評估數(shù)值實驗和時間數(shù)值實驗的所有程序都是以SCILAB代碼撰寫的。數(shù)值實驗的硬件平臺是Intel(R) Core(TM) i5-2450M,CPU主頻為2.50 GHz，內(nèi)存為4.00 GB和飛騰處理器FT1500A，主頻為1.5 GHz，內(nèi)存為32 GB；所用的操作系統(tǒng)是UbantuKylin 14.04；開源軟件版本為SCILABv 5.5.2。

3.1 冪指數(shù)基多項式函數(shù)值的高精度可靠試驗

多項式的標(biāo)準(zhǔn)格式為：

p(x)=a0+a1x+…+anxn

工程和數(shù)學(xué)的許多領(lǐng)域需要求多項式在某一點的函數(shù)值或?qū)?shù)值，最常用的是Horner算法[16]。Horner 算法能夠同時計算多項式函數(shù)的任意階導(dǎo)數(shù)，Müller[17]、Olver[18]等學(xué)者對Horner 算法進(jìn)行了深入研究，并給出了數(shù)值誤差的理論分析。經(jīng)過研究發(fā)現(xiàn)，Horner算法具有較好的穩(wěn)定性，并且在長時程、大規(guī)模甚至是病態(tài)的數(shù)據(jù)運算中，Horner算法都具有比較明顯的優(yōu)勢。

Horner算法解多項式的過程如算法3所示。

算法3Horner

輸入：p,x。

輸出：q0。

步驟1qn=an；

步驟2 fori=n-1:-1:0then

步驟3qi(x) =xqi+1(x)+ai；

步驟4q0=q0(x)。

基于無誤差變換，同樣寫出double-double格式下的Horner算法，即Horner_DD(p,x),如算法4所示。

算法4Horner_DD

輸入：p,x。

輸出：res。

步驟1x=DD(x)；

步驟2y=DD(0)；

步驟3 fori=n:-1:1then

步驟4y=x*y+a(i);

步驟5res=hi(y)。

實驗使用Horner算法對多項式p(x)=(x-2)9進(jìn)行展開，并計算在其重根x=2附近區(qū)域的值。在區(qū)間[1.92,2.08]中均勻取樣8 000個點，計算得到的結(jié)果如圖1所示，分別使用double和double-double數(shù)據(jù)格式計算得出。從圖1中可以看到，Horner_DD算法與Horner算法相比，繪出的曲線更加光滑，得到了更好的計算結(jié)果。

Figure 1 Horner rules in double and double-double format圖1 應(yīng)用double和double-double數(shù)據(jù)格式的Horner算法繪制的多項式展開式的圖形

經(jīng)統(tǒng)計，Horner算法的運行時間time_d=0.2160，Horner_DD算法的運行時間time_dd=10.78939，兩者運行時間比為time_dd/time_d=34.87。Horner算法共執(zhí)行2n次浮點運算，Horner_DD算法共執(zhí)行34n次浮點運算。后者運算量是前者的17倍，計算時間是前者的34.87倍。效率變慢的原因是更多的函數(shù)調(diào)用和加載。

3.2 Bernstein基多項式函數(shù)的高精度可靠實驗

計算機(jī)輔助幾何設(shè)計中，單變量多項式一般表示為Bernstein基形式,曲線p(t)的定義為：

求解曲線p(t)的算法為De Casteljau算法，簡稱DC算法,如算法5所示。

算法5DC

輸入：p,t。

輸出：res。

步驟1b(i,1)=b(i),i=1,…,n；

步驟2 forj=1:1:nthen

步驟3fori=0:1:n-jthen

步驟4b(i+1,j+1)=(1-t)*b(i+1,j)+t*b(i+2,j)；

步驟5res=b(1,n+1)。

基于無誤差變換，同樣寫出double-double格式下的DC算法，表示為DC_DD(p,t),如算法6所示。

算法6DC_DD

輸入：p,t。

輸出：res。

步驟1t=DD(t)；

步驟2b(i)(1)=list(p(i),0),i=1,…,n；

步驟3 forj=1:1:n-1then

步驟4fori=0:1:n-1-jthen

步驟5b(i+1)(j+1)=(1-t)*b(i+1)(j)+t*b(i+2)(j)；

步驟6temp=b(1)(n)；

步驟7res=hi(temp)。

Figure 2 DC rules in double and double-double format圖2 應(yīng)用double和double-double數(shù)據(jù)格式的DC算法繪制多項式展開式的圖形

經(jīng)統(tǒng)計，DC算法的運行時間time_d=1.460，DC_DD算法的運行時間time_dd=54.76,兩者運行時間比為time_dd/time_d=44.73。DC算法的浮點計算量為1.5n2+1.5n+1，DC_DD算法的浮點計算量為33n2+33n+2，后者的運行時間是前者的44.73倍，浮點計算量是前者22倍。效率變慢的原因同樣是更多函數(shù)加載和調(diào)用。

3.3 Chebyshev基多項式函數(shù)的高精度可靠實驗

Chebyshev多項式是正交多項式的一種，由于其在逼近論中的優(yōu)異特性而得到廣泛應(yīng)用。計算Chebyshev多項式在某一點的函數(shù)值采用迭代的Clenshaw算法[19]。Chebyshev多項式為定義在[-1,1]上的經(jīng)典正交多項式。使用Chebyshev基對多項式r(x)展開為：

其中，bi為展開系數(shù)，Ai(x)為Chebyshev多項式：

A0(x)=1

A1(x)=x

An+1(x)=2xAn(x)-An-1(x)

Clenshaw算法描述如算法7所示。

算法7Clenshaw

輸入：r,x。

輸出：res。

步驟1an+2=an+1=0；

步驟2 forj=n:-1:1then

步驟3aj=2xaj+1-aj+2+Chebj；

步驟4res=xa1-a2+b0。

基于無誤差變換，同樣寫出double-double格式下的Clenshaw算法，表示為Clenshaw_DD(r,x),如算法8所示。

算法8Clenshaw_DD

輸入：r,x。

輸出：res。

步驟1x=DD(x)；

步驟2DD(a(i))；

步驟3 forj=n:-1:2then

步驟4a(j)=2*x*a(j+1)-a(j+2)+Cheb(j);

步驟5temp=x*a(2)-a(3)+Cheb(1)；

步驟6res=hi(temp)。

實驗計算多項式r(x)=(x-0.75)7(x-1)10的Chebyshev基展開。本文選取一組17項展開的病態(tài)數(shù)據(jù)bi，然后在[0.68,1.15]上取樣8 000個均勻點，分別用Clenshaw算法和Clenshaw_DD算法繪出曲線，結(jié)果如圖3所示。

從圖3可以看出，同Clenshaw算法相比，Clenshaw_DD算法給出了較為光滑的曲線，有更好的計算結(jié)果。

Figure 3 Clenshaw rules in double and double-double format圖3 應(yīng)用double和double-double數(shù)據(jù)格式Clenshaw算法繪制多項式展開式的圖形

經(jīng)統(tǒng)計，Clenshaw算法所用平均時間time_d=0.876，Clenshaw_DD算法所用平均時間time_dd=30.28。兩者運行平均時間之比time_dd/time_d=43.69。Clenshaw算法的浮點計算量為3n+4，Clenshaw_DD算法的浮點計算量為52n+53。采用double-double格式，在計算量提高14倍的情況下，運行時間增加了43.49倍，造成效率低效的原因是更多的函數(shù)加載和調(diào)用。

從前面的實驗結(jié)果可以看到，使用double-double數(shù)據(jù)格式設(shè)計多精度算法庫，雖然增加了浮點運算量，但實現(xiàn)了對舍入誤差的有效控制，數(shù)值計算結(jié)果的可靠性得到了保障。

3.4 比較與分析

首先驗證3個數(shù)值驗證實驗成果的正確性和有效性。在CPU主頻為2.50 GH、內(nèi)存為4.00 GB、操作系統(tǒng)為Ubuntu Kylin 14.04的Intel處理器上統(tǒng)計CPU的運行時間，得到表5。

Table 5 Run time of algorithms on Intel processor

從表5可以看出，時間比time_dd/time_d約為30多倍。double-double格式更多地記錄了浮點數(shù)的舍入誤差，計算量更大，調(diào)用函數(shù)更多，所以花費的時間也比double格式的多。

在Intel處理器上構(gòu)建好函數(shù)庫后，將其移植到飛騰處理器平臺上，統(tǒng)計CPU運行時間，得到表6。

Table 6 Run time of algorithms on FT processor

從表6可以看出，時間比time_dd/time_d約為30多倍。與在Intel處理器上的結(jié)果相近，double-double格式更多地記錄了浮點數(shù)的舍入誤差，計算量更大，調(diào)用函數(shù)更多，所以花費的時間也比double格式的多。

計算出飛騰處理器上double-double與double時間的比值，與Intel處理器的數(shù)據(jù)進(jìn)行比較，得到圖4。

Figure 4 Run time ratio of algorithmsin Intel processor and FT processor圖4 算法在Intel處理器和飛騰處理器上運行的時間比time_dd/time_d

2種處理器上，比值time_dd/time_d大致相同，這表明double-double格式下構(gòu)建的高精度算法庫具有較好的可移植性。

4 結(jié)束語

隨著科學(xué)工程計算規(guī)模的增加、計算時程的增長，數(shù)值結(jié)果的精度可靠性越來越受到研究者的重視。對于國防安全等需求，如何在國產(chǎn)處理器上實現(xiàn)自主可控的軟件至關(guān)重要。本文面向開源軟件平臺SCILAB，設(shè)計了多精度數(shù)值算法庫。數(shù)值實驗結(jié)果表明，雖然多精度算法增加了一定的浮點運算量，但計算結(jié)果的可靠性得到了進(jìn)一步的保證，有效控制了計算機(jī)舍入誤差的累積。與此同時，通過對國產(chǎn)自主可控處理器上的結(jié)果和Intel處理器上的結(jié)果對比，表明本文設(shè)計的多精度算法庫具有較好的可移植性。本文的工作可以用于提升大規(guī)?？茖W(xué)工程計算的穩(wěn)定性和可靠性，以及為在國產(chǎn)自主可控處理器上實現(xiàn)自主可控軟件提供了支持。未來準(zhǔn)備在前期工作的基礎(chǔ)上，繼續(xù)豐富該算法庫的功能，為應(yīng)用提供更多支持。