基于SCILAB的多精度算法研究與實現 *

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

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

1 引言

高性能處理器被廣泛應用于信息化聯合作戰、航空航天、模擬核試驗、氣象預報和搶險救災等領域。近年來，我國自主研發的處理器發展迅速，這對于我國的國家安全具有重要意義，有利于完善我國經濟信息系統以及國防領域的安全和保密問題。隨著計算機技術的不斷發展，高端計算已經開始從單一追求高性能轉變為致力于綜合高效能的提升，以解決當前高性能計算領域所面臨的實用性能和可編程性問題。浮點數運算對于計算機性能影響顯著，在浮點運算過程中存在舍入誤差，計算規模較小時，這些誤差不會造成明顯影響。現代科學工程計算一般涉及到大規模計算，在這種大規模、長時程的數值計算中，由于舍入誤差的累積效應存在，可能導致數值計算結果不準確、不可靠，甚至得到完全錯誤的結果。高精度算法庫的設計和實現，是解決大規模數值計算的可靠性和穩定性，以及提升國產處理器性能的有效途徑之一。

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

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

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

2 多精度算法研究與實現

2.1 無誤差變換

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

(a°b)=x+y

(1)

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

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

算法1FastTwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a⊕b；

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

(2)兩個浮點數乘法的無誤差變換[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個浮點數拆成1對浮點數。

2.2 double-double算術

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

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

x=xh+xl

其中，xh和xl沒有交疊，且xh和xl是符合IEEE-754標準的雙精度double格式數，那么x就是double-double格式數。

2.3 運算符重載

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

二元操作數變換規則：

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

一元操作數變換規則：

%〈operand_type〉_〈op_code〉

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

本文將對表3中所列出的運算符進行重載。

根據運算符重載的規則，以加法重載的實現為

Table 1 Cross reference list to operand_type

Table 2 Cross reference list to op_code

Table 3 Overloading operational character

例，對運算符重載的過程進行展示。把浮點數a，b分割成高位H和低位L，分別記為ah，al，bh，bl，然后實現double-double的加法。在此基礎上使用重載命令對加法進行重載。基于double-double數據格式的加法重載過程為：

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 函數重載與構建

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

Table 4 Overloading functions

以equal函數為例，本文基于double-double格式對其進行重載，重載具體實現過程為：

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格式重新定義的函數，使用deff進行重載操作。

3 數值實驗與分析

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

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

多項式的標準格式為：

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

工程和數學的許多領域需要求多項式在某一點的函數值或導數值，最常用的是Horner算法[16]。Horner 算法能夠同時計算多項式函數的任意階導數，Müller[17]、Olver[18]等學者對Horner 算法進行了深入研究，并給出了數值誤差的理論分析。經過研究發現，Horner算法具有較好的穩定性，并且在長時程、大規模甚至是病態的數據運算中，Horner算法都具有比較明顯的優勢。

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進行展開，并計算在其重根x=2附近區域的值。在區間[1.92,2.08]中均勻取樣8 000個點，計算得到的結果如圖1所示，分別使用double和double-double數據格式計算得出。從圖1中可以看到，Horner_DD算法與Horner算法相比，繪出的曲線更加光滑，得到了更好的計算結果。

Figure 1 Horner rules in double and double-double format圖1 應用double和double-double數據格式的Horner算法繪制的多項式展開式的圖形

經統計，Horner算法的運行時間time_d=0.2160，Horner_DD算法的運行時間time_dd=10.78939，兩者運行時間比為time_dd/time_d=34.87。Horner算法共執行2n次浮點運算，Horner_DD算法共執行34n次浮點運算。后者運算量是前者的17倍，計算時間是前者的34.87倍。效率變慢的原因是更多的函數調用和加載。

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

計算機輔助幾何設計中，單變量多項式一般表示為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 應用double和double-double數據格式的DC算法繪制多項式展開式的圖形

經統計，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倍。效率變慢的原因同樣是更多函數加載和調用。

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

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

其中，bi為展開系數，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項展開的病態數據bi，然后在[0.68,1.15]上取樣8 000個均勻點，分別用Clenshaw算法和Clenshaw_DD算法繪出曲線，結果如圖3所示。

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

Figure 3 Clenshaw rules in double and double-double format圖3 應用double和double-double數據格式Clenshaw算法繪制多項式展開式的圖形

經統計，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倍，造成效率低效的原因是更多的函數加載和調用。

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

3.4 比較與分析

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

Table 5 Run time of algorithms on Intel processor

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

在Intel處理器上構建好函數庫后，將其移植到飛騰處理器平臺上，統計CPU運行時間，得到表6。

Table 6 Run time of algorithms on FT processor

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

計算出飛騰處理器上double-double與double時間的比值，與Intel處理器的數據進行比較，得到圖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格式下構建的高精度算法庫具有較好的可移植性。

4 結束語

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