趙永浩,賈海鵬,張云泉,張思佳

(1.大連海洋大學(xué)信息工程學(xué)院，遼寧 大連 116023;2.中國科學(xué)院計算技術(shù)研究所計算機(jī)體系結(jié)構(gòu)國家重點(diǎn)實(shí)驗(yàn)室，北京100190)

1 引言

平方根函數(shù)是基本的數(shù)學(xué)運(yùn)算，也是處理器軟件生態(tài)的重要組成部分。平方根函數(shù)的應(yīng)用極其廣泛[1]，可應(yīng)用于圖形處理器[2]、神經(jīng)網(wǎng)絡(luò)、計算機(jī)圖形學(xué)和積分計算等科學(xué)領(lǐng)域。隨著處理器架構(gòu)的不斷發(fā)展和完善，為了更好地支持和完善計算機(jī)的并行性,SIMD(Single Instruction Multiple Data)技術(shù)得到了快速的發(fā)展，已經(jīng)被廣泛使用到各種處理器上，并成為其中重要的技術(shù),比如ARM架構(gòu)的NEON指令集[3]，MIPS的MSA指令集和Intel x86架構(gòu)中的MMX/SSE/AVX指令集[4]。NEON技術(shù)是對ARM Cortex-A和Cortex-R系列處理器中SIMD的擴(kuò)展。NEON寄存器存儲的是同一數(shù)據(jù)類型的向量，相比其他指令其優(yōu)點(diǎn)是NEON指令可以對4個32 bit和2個64 bit的元素同時進(jìn)行計算，并且還支持int、float、double等多種數(shù)據(jù)類型。該技術(shù)通過同時計算多個元素，加快計算速度，來達(dá)到提升整體性能的目的。

研究基于SIMD的平方根函數(shù)的高性能實(shí)現(xiàn)與優(yōu)化具有十分重要的研究價值和實(shí)際意義。針對ARM架構(gòu)處理器在市場供不應(yīng)求的現(xiàn)狀，SIMD向量化技術(shù)也在不斷發(fā)展，所以完善ARM SIMD生態(tài)是十分重要的工作。目前基于SIMD NEON技術(shù)的高性能函數(shù)實(shí)現(xiàn)與優(yōu)化的研究比較少，可以查閱到的資料有對處理器的擴(kuò)展函數(shù)庫的實(shí)現(xiàn)與優(yōu)化，另一個是針對SIMD功能的函數(shù)優(yōu)化算法，但這2項研究的時間較早，技術(shù)還不算完善。所以，利用SIMD NEON技術(shù)針對ARM架構(gòu)進(jìn)行平方根函數(shù)高性能實(shí)現(xiàn)具有重大意義。平方根函數(shù)高性能實(shí)現(xiàn)與優(yōu)化[5]的基本思想是，在保證精度的條件下，盡可能地提升算法的性能。所以，設(shè)計高精度、高性能的平方根函數(shù)算法是本文的研究重點(diǎn)。就目前來說，處理器已經(jīng)提供了sqrt指令，可實(shí)現(xiàn)平方根的SIMD求解，但其性能較低。在測試數(shù)據(jù)長度相同的情況下，直接調(diào)用SIMD求平方根的指令，計算時間需要0.214 570 ms，而本文的計算方法只需要0.069 750 ms，性能可以提高將近3倍。平方根函數(shù)計算方法是基于數(shù)值分析[6]領(lǐng)域的方法，其中用到的算法是平方根倒數(shù)和泰勒展開[7]，使用2個算法的目的是使計算在滿足計算精度的同時盡量提升計算時間性能。本文實(shí)驗(yàn)中，將與本文算法與libm算法庫[8]進(jìn)行性能比較，并和處理器本身提供的sqrt指令進(jìn)行比較。本文的主要貢獻(xiàn)如下所示：

(1)依據(jù)計算平方根函數(shù)的方法，得出平方根函數(shù)高性能實(shí)現(xiàn)的思想和原理；

(2)總結(jié)出平方根函數(shù)高性能實(shí)現(xiàn)和優(yōu)化的具體方法；

(3)實(shí)現(xiàn)了具有高性能平方根函數(shù)的算法，并與libm算法庫進(jìn)行了性能對比。

2 相關(guān)工作

本文主要基于數(shù)值分析的相關(guān)原理利用SIMD對平方根函數(shù)的高性能實(shí)現(xiàn)進(jìn)行優(yōu)化；現(xiàn)有的一些函數(shù)庫提供了計算平方根函數(shù)的算法，如ARM_M算法庫、ARM Performance Library(ARMPL)和libm算法庫。本文最后要與libm算法庫和處理器提供的sqrt指令在計算性能和精度2個方面進(jìn)行對比。

2.1 libm

Newlib是一個面向嵌入式系統(tǒng)的C運(yùn)行庫，對于與GNU兼容的嵌入式C運(yùn)行庫，Newlib并不是唯一的選擇，但Newlib具有獨(dú)特的體系結(jié)構(gòu)，使得它能夠非常好地滿足嵌入式系統(tǒng)的要求，此外，Newlib具有很強(qiáng)的可移植性和可重入性，且功能完備，已廣泛應(yīng)用于各種嵌入式系統(tǒng)中。

libm算法庫是C運(yùn)行庫Newlib的一部分，也是Newlib中一個完整的IEEE數(shù)學(xué)庫，提供了用于嵌入式系統(tǒng)接收float和double 2種類型參數(shù)的函數(shù)。

在基于 SPARC 的系統(tǒng)上，libm 中的初等函數(shù)是使用表驅(qū)動算法和多項式有理數(shù)近似算法的組合實(shí)現(xiàn)的，可以實(shí)現(xiàn)更高的性能或精確度。

在基于x86的系統(tǒng)上，libm的某些初等函數(shù)是使用x86指令集中提供的初等函數(shù)內(nèi)核指令實(shí)現(xiàn)的；其他函數(shù)是用基于 SPARC 的系統(tǒng)上的相同表驅(qū)動算法或多項式/有理數(shù)近似算法實(shí)現(xiàn)的。

libm 中常見單精度初等函數(shù)的表驅(qū)動算法和多項式/有理數(shù)近似算法都將準(zhǔn)確的計算結(jié)果傳送到最后位置單元內(nèi)，通過運(yùn)算法則直接分析這些誤差界限。

2.2 ARMPL

ARMPL是ARM公司為滿足科學(xué)計算和高性能計算社區(qū)應(yīng)用開發(fā)需求所推出的，針對ARM V8-A系列處理器精準(zhǔn)化的商業(yè)數(shù)學(xué)函數(shù)庫。ARMPL通過充分利用芯片合作伙伴為ARM V8-A架構(gòu)系統(tǒng)單芯片設(shè)計的創(chuàng)新特性和功能，確?；贏RM V8-A架構(gòu)的高性能服務(wù)器和系統(tǒng)到達(dá)極致，在科學(xué)計算應(yīng)用程序使用中，可以確保獲得的結(jié)果更加準(zhǔn)確。

2.3 ARM V8 架構(gòu)

ARM 是一種負(fù)載性的存儲體系結(jié)構(gòu)，是RISC 處理器的典型之一。ARM V8兼容了ARM V7，是一個支持64 位指令集的ARM 處理器架構(gòu)，引入64 位體系結(jié)構(gòu)之后，還保持了與現(xiàn)有32 位系統(tǒng)結(jié)構(gòu)的兼容性。

ARM V8 提供了31個64 bit 通用寄存器和32 個128 bit 浮點(diǎn)寄存器，浮點(diǎn)寄存器在執(zhí)行指令時一條指令可以操作多個操作數(shù)，可以提高指令的執(zhí)行效率。ARM V8架構(gòu)提供了sqrt指令，可以直接對浮點(diǎn)數(shù)進(jìn)行開平方計算。本文將與sqrt指令在性能方面進(jìn)行對比。

2.4 向量化SIMD

SIMD用一個控制器來操作多個處理單元，時間上的并行性是通過對數(shù)據(jù)中的每個元素進(jìn)行相同的計算來實(shí)現(xiàn)的。目前該技術(shù)已經(jīng)在多個平臺使用，而且每個平臺都擁有不同的指令集。編譯器把代碼翻譯成SIMD指令的操作叫做SIMD并行化或者SIMD向量化。在執(zhí)行SIMD指令時，通過一條指令可以同時對多個元素進(jìn)行計算。例如，在一個128 bit的向量寄存器下，每次編譯時，可以對4個float類型數(shù)據(jù)或2個double類型數(shù)據(jù)進(jìn)行計算。理論上，在性能方面，SIMD 指令是普通運(yùn)算的4倍。

2.5 平方根算法綜述

目前計算平方根的方法有以下幾種：Newton-Raphson(牛頓一拉斐森)算法[9]、Goldschmitdt算法、SRT-Redundant算法[10]和Non-Redundant算法[11]。

Newton-Raphson算法[9]已經(jīng)廣泛應(yīng)用于許多工程。該算法的實(shí)質(zhì)是利用泰勒展開定理把非線性方程展開為線性方程，在迭代過程中需要使用乘法操作來計算除數(shù)的倒數(shù)，這會導(dǎo)致延時變大，此外，迭代的周期數(shù)與初始值有關(guān)，所以沒有固定的周期。

Goldschmitdt算法是Newton-Raphson迭代算法的改進(jìn)，通過增加子操作之間的并行性來減小延時。

上述2種算法的缺點(diǎn)是運(yùn)算速度慢，在精度上無法滿足IEEE浮點(diǎn)數(shù)標(biāo)準(zhǔn)規(guī)定的精度要求。

Non-Redundant算法[11]與SRT-Redundant算法類似，但該算法是用補(bǔ)碼的形式來表示平方根，其運(yùn)算精度低。

3 研究背景

本文平方根函數(shù)主要使用的是平方根倒數(shù)算法和泰勒展開[12]。

3.1 平方根倒數(shù)速算法

IEEE 754標(biāo)準(zhǔn)浮點(diǎn)數(shù)在內(nèi)存中的存儲方式如圖1所示。

Figure 1 IEEE 754 float point format圖1 IEEE 754浮點(diǎn)數(shù)格式

IEEE 754標(biāo)準(zhǔn)中,一個規(guī)格化浮點(diǎn)數(shù)的真值可以表示為式(1)所示：

x=(-1)s*(1.M)*2e

(1)

其中,s代表符號位，e(exponent)代表指數(shù)，單精度浮點(diǎn)數(shù)e=E-127，雙精度浮點(diǎn)數(shù)e=E-1023,1.M代表尾數(shù)域，

平方根倒數(shù)速算法(Fast Inverse Square Root)[13]是適用于快速計算積的算術(shù)平方根的倒數(shù)的一種算法，公式如式(2)所示：

(2)

其中，x是符合IEEE 754標(biāo)準(zhǔn)格式的32位浮點(diǎn)數(shù)(該算法是為IEEE 754標(biāo)準(zhǔn)中的浮點(diǎn)數(shù)存儲規(guī)則設(shè)計的)。該算法的優(yōu)勢在于減少了求平方根倒數(shù)時浮點(diǎn)運(yùn)算操作帶來的巨大運(yùn)算消耗，大大提高了運(yùn)算速度，比調(diào)用sqrt指令快了約20倍。但是，該算法所得到的結(jié)果精度不夠，所以本文將采用泰勒展開公式對其精度進(jìn)行補(bǔ)償。

3.2 泰勒級數(shù)

如果函數(shù)f(x)在點(diǎn)x0處可以求導(dǎo)，則有:

f(x)=f(x0)+f′(x0)(x-x0)+o(x-x0)

(3)

當(dāng)|x-x0|的結(jié)果很小時，則有：

f(x)≈f(x0)+f′(x0)(x-x0)

(4)

等式右邊是等式左邊的近似值，但右邊的式子明顯精度不高，與等式左邊相比誤差過大。所以，可將上式推廣至高次多項式：

f(x)=f(x0)+f′(x0)(x-x0)+

(5)

等式右邊的多項式次數(shù)越高，對應(yīng)結(jié)果的精度也越高。但是，若要達(dá)到可以滿足條件的精度，有時相應(yīng)的多項式次數(shù)會過高，造成計算時間增加到無法接受的地步。

(6)

4 平方根函數(shù)的算法實(shí)現(xiàn)

本文研究了4種類型的平方根函數(shù)的計算方法，包括單精度浮點(diǎn)數(shù)、雙精度浮點(diǎn)數(shù)、單精度復(fù)數(shù)和雙精度復(fù)數(shù)，每一種類型函數(shù)的實(shí)現(xiàn)都是基于C語言的代碼進(jìn)行SIMD并行優(yōu)化[15]。下面主要介紹計算每種數(shù)據(jù)類型平方根函數(shù)的算法的高性能實(shí)現(xiàn)的原理及其計算方法。

4.1 單精度浮點(diǎn)數(shù)的算法實(shí)現(xiàn)

就目前來說，處理器已經(jīng)提供了計算sqrt的指令實(shí)現(xiàn)平方根的SIMD求解，但是在性能方面比較差。在測試數(shù)據(jù)長度相同的情況下，直接調(diào)用SIMD求平方根的指令，計算時間需要0.214 570 ms，而本文方法只需要0.069 750 ms，性能可以提高將近3倍。

(7)

(8)

(9)

將(x/x0)-0.5根據(jù)泰勒公式展開,得：

(10)

最終表達(dá)式為：

(11)

在計算過程中，用平方根倒數(shù)的算法得出的結(jié)果只能精確到前8 bit，所以本文用泰勒展開對精度進(jìn)行補(bǔ)償,來確保計算精度。具體計算過程如算法1所示：

算法1單精度浮點(diǎn)數(shù)的計算

Input:x0。//約等于x的數(shù)

Output:y。

1．a(chǎn)←rsqrt(x)；

2．b←x*a；

3．y5←0.5*a；

4．y9←0.5 *(1-ab)；

5．y11←b*(1+0.5*(1-ab))；

6．y5←0.5*a(1+0.5*(1-ab))；

7．y←y11+(x-y112)*y5。

通過該計算過程計算得出的結(jié)果，在性能和精度方面相比直接調(diào)用sqrt指令有了很大的提升，性能方面可以提高大約7倍；精度方面，平方根倒數(shù)算法可以精確到前8 bit，此后用泰勒展開的方法對精度進(jìn)行補(bǔ)償，2種方法結(jié)合在一起，既保證了性能，也達(dá)到了對精度的要求。

4.2 雙精度浮點(diǎn)數(shù)的算法實(shí)現(xiàn)

(12)

具體包括如下步驟：

步驟2對尾數(shù)域1.M進(jìn)行計算，如果直接對1.M開根號計算，不僅計算的時間較長，而且計算結(jié)果會有很大誤差。所以，本文將1.M拆分成2個單獨(dú)的部分分別進(jìn)行計算。令1.M=x，然后將x表示為x0和x1，對x0進(jìn)行計算用的是開根號取倒數(shù)的方法，此方法只針對單精度浮點(diǎn)數(shù)；對x1進(jìn)行計算用的根據(jù)泰勒展開公式的方法。具體算法將在后面介紹。

(13)

(14)

最終的計算公式為：

(15)

需要使用聯(lián)合體(union)來把輸入重新解釋為浮點(diǎn)數(shù)、整型數(shù)和無符號整型數(shù)，并且通過聯(lián)合體可以使用邏輯運(yùn)算求得浮點(diǎn)數(shù)的指數(shù)位和尾數(shù)位。

使用union將雙精度浮點(diǎn)數(shù)x重新解釋為一個無符號的整數(shù)，通過邏輯運(yùn)算取出指數(shù)位，以判斷浮點(diǎn)數(shù)x指數(shù)的奇偶性。下面介紹對指數(shù)位的計算過程：

(1)首先，為了將11位的指數(shù)位取出，利用union將64位浮點(diǎn)數(shù)重新解釋為64位無符號的整數(shù)。

(2)用0x7ff000000000000與重新解釋為64位無符號的整數(shù)進(jìn)行“與”運(yùn)算，通過運(yùn)算可以得到指數(shù)位e，其他位置0。

(3)接下來判斷指數(shù)的奇偶性，將取出后的指數(shù)位按位取反,將～e和0x0010000000000000進(jìn)行一次“與”運(yùn)算。

(4)將(3)的計算結(jié)果加上0xbfd0000000000000,然后用該結(jié)果再減去(2)取得的指數(shù)位。

通過以上運(yùn)算，得出e的表達(dá)式為：

(1)當(dāng)e% 2=0時：

(16)

(2)當(dāng)e% 2=1時：

(17)

接下來介紹對x0進(jìn)行的運(yùn)算。計算x0是按照單精度浮點(diǎn)數(shù)計算方式，在計算指數(shù)部分時，先對取出來的指數(shù)位進(jìn)行奇偶性判斷，然后再與x0整合在一起，形成一個新的單精度浮點(diǎn)數(shù)。計算x0的具體過程如下所示：

(1)設(shè)e是浮點(diǎn)數(shù)的指數(shù)位，將e與0x07f0000000000000相加,得到相加結(jié)果后再整體往右移動29位。

(2)用0x001fffffffffffff和重新解釋為64位無符號的整數(shù)進(jìn)行“與”運(yùn)算，得到尾數(shù)位，其他位置0。然后將尾數(shù)位整體往右移動29位，這樣雙精度浮點(diǎn)數(shù)的尾數(shù)位上前29位為0，后23位為移動前尾數(shù)位的高23位，將其當(dāng)作新浮點(diǎn)數(shù)的尾數(shù)位。

(3)將(1)和(2)的結(jié)果進(jìn)行“或”運(yùn)算后，組成一個新的32位單精度浮點(diǎn)數(shù)，然后就可以用平方根倒數(shù)的算法對其進(jìn)行計算。

(4)平方根倒數(shù)算法得出的結(jié)果為float類型，需要將其轉(zhuǎn)換成double類型，可使用SIMD NEON提供的verinterpretq_f64_u64轉(zhuǎn)換指令完成。這條指令在轉(zhuǎn)換后沒有改變值，只改變數(shù)據(jù)類型。具體計算過程如算法2所示。

算法2x0計算過程

Input:x。

Output:y0。

1.DECLAREsrcAS DOUBLE_U；

2.DECLAREsAS FLOAT_U；

3.DECLAREtempAS FLOAT；

4.DECLAREe,src.uUNSIGNED INT；

5.src.f←x；

6.e←src.u&x_index；

7.e←0x001000000000000 & (～e)；

8.e←e+0x07f0000000000000；

9.e←e>> 29；

10.u←src.u&x_mantissa；

11.u←u>> 29；

12.s.u←u|e；

13.temp←rsqrt(s.f)；

14.y0←temp。

在計算x1時使用的是泰勒公式展開的方法，具體計算過程如下所示:

(1)在取出雙精度浮點(diǎn)數(shù)的尾數(shù)位的運(yùn)算中，使用0x001fffffffffffff和64位無符號整數(shù)進(jìn)行“與”運(yùn)算,得到的結(jié)果會保留指數(shù)位中最后一位。在下面的計算中采用位移的方法將指數(shù)位中末位的1去掉。

x1=1-x*y2

(18)

(19)

具體計算過程如算法3所示。

算法3x1計算過程

Input:x,y0。

Output:dst。

1.DECLAREb,dAS DOUBLE；

2.DECLAREe,u,EUNSIGNED INT；

3.DECLAREsrc,sAS DOUBLE _U；

4.src.f←x；

5.e←src.u&x_index；

6.e←0x001000000000000 & (～e)；

7.E←e+0xbfd0000000000000；

8.u←E-e；

9.s.u←u>> 1；

10.y←s.f*y0；

11.b←y*y；

12.x1←1-b*x；

13.d←x1*(x1*((x1*A5)+A4)+A3)+A2；

以上就是計算雙精度浮點(diǎn)數(shù)的方法和過程，在保證雙精度浮點(diǎn)數(shù)所能表示的數(shù)據(jù)精度和計算速度的同時，完成浮點(diǎn)數(shù)的開平方運(yùn)算。

4.3 復(fù)數(shù)的平方根算法實(shí)現(xiàn)

在對復(fù)數(shù)進(jìn)行平方根計算的時候，本文使用的是將復(fù)數(shù)的實(shí)部和虛部分開計算的方法。

(20)

(21)

通過2個公式可以分別求出開根號之后的實(shí)部和虛部結(jié)果。具體計算過程如算法4所示。

算法4復(fù)數(shù)平方根的計算

Input：z=a+bi。

Output：w=c+di。

1．a(chǎn)=p_src.re；

2．b=p_src.im；

3．R=sqrt(a2+b2)；

4．c=sqrt((a+R)*0.5)；

5．temp.d=b*0.5；

6．d=temp.d/c。

5 SIMD優(yōu)化

根據(jù)第4節(jié)函數(shù)算法的實(shí)現(xiàn)，在執(zhí)行平方根函數(shù)的計算時，可以利用ARM NEON INTRINSIC指令對SIMD代碼進(jìn)行優(yōu)化，以此加快對平方根函數(shù)的計算。此外，由于ARM NEON INTRINSIC提供的SIMD乘加指令的速度很快，因此可以在對程序進(jìn)行精度補(bǔ)償時使用泰勒展開公式，從而使用多個SIMD乘加指令加速運(yùn)算。

在華為的鯤鵬920平臺下，每個向量寄存器的長度為128 bit，所以可以同時計算4個float類型數(shù)據(jù)或者2個double類型數(shù)據(jù)。在ARM NEON INTRINSIC提供的指令中，只有l(wèi)oad和store是寄存器與存儲器之間的操作，其它指令都是介于寄存器與寄存器之間的操作。所以，可以把所有的計算放到寄存器中完成。

此外，在計算單精度浮點(diǎn)數(shù)時直接用了SIMD NEON的rsqrt指令，指令的功能是計算輸入x的平方根再取倒數(shù)的值。此外，在計算雙精度浮點(diǎn)數(shù)的過程中，采用了循環(huán)展開的方法，因?yàn)橐粋€double類型數(shù)據(jù)是64 bit，而一個寄存器的長度為128 bit，所以用循環(huán)展開的方式，一次可以計算4個雙精度浮點(diǎn)數(shù)，然后通過使用2個store把4個結(jié)果分別放到2個存儲器中，這樣可以加快計算速度。

6 性能評估

6.1 測試環(huán)境

(1)硬件環(huán)境。

本文實(shí)驗(yàn)采用鯤鵬920作為測試平臺，鯤鵬920使用的是ARM架構(gòu)，主頻為2.6 GHz，測試環(huán)境的配置如表1所示。

Table 1 Hardware environment of experiment

(2)軟件環(huán)境。

本文軟件環(huán)境使用ARM V8架構(gòu)提供的sqrt指令和C語言標(biāo)準(zhǔn)算法庫作為計算速度的對比對象。對于單精度浮點(diǎn)數(shù)和雙精度浮點(diǎn)數(shù)的測試樣例輸入，輸出的平均精度需小于1E-6,最大誤差小于1E-15才可以通過樣例測試。

6.2 性能分析

如表2 所示，Open-v 是本文提出的算法實(shí)現(xiàn)，libm是C 語言標(biāo)準(zhǔn)算法庫，sqrt 是ARM V8架構(gòu)提供的求平方根函數(shù)的指令。表2中性能指標(biāo)測試的輸入向量長度為216。性能指標(biāo)的單位是ns/item，即每個輸入值計算所需要的平均時間。

Table 2 Performance comparison

由表2可以看到,本文工作相對于libm算法庫在性能方面具有很大的提升，比libm算法庫快了約7倍，比調(diào)用sqrt指令快了約4倍。并且可以看出算法對單精度函數(shù)的優(yōu)化效果相較于雙精度函數(shù)更明顯。原因是計算float類型數(shù)據(jù)時可以直接調(diào)用指令來計算，縮短了計算時間。而在計算double類型數(shù)據(jù)時，算法比較復(fù)雜，首先，需要將尾數(shù)中的前23位取出來，當(dāng)作float類型數(shù)據(jù)來計算；然后，再將計算結(jié)果轉(zhuǎn)換成double類型數(shù)據(jù)；最后，需要利用泰勒展開公式對其展開多次，以達(dá)到對計算結(jié)果的精度要求。這是在算法上造成的性能差異的主要原因。另外的原因是在內(nèi)存中的表示形式，float類型數(shù)據(jù)占4 B(32 bit)內(nèi)存空間，double類型數(shù)據(jù)占8 B(64 bit)內(nèi)存空間。在一個時鐘周期內(nèi)，CPU發(fā)射一個計算指令可以計算2個float類型數(shù)據(jù)，而double類型數(shù)據(jù)只能計算1個。綜上所述，計算double類型數(shù)據(jù)要比計算float類型數(shù)據(jù)在性能上差距更大。

在性能指標(biāo)測試的輸入向量長度為65 536時，直接調(diào)用ARM V8架構(gòu)提供的求平方根函數(shù)的sqrt指令計算時間為3 383 ns，而直接調(diào)用rsqrt指令計算時間為185 ns，可以看出速度比sqrt快約20倍，但是調(diào)用rsqrt指令計算，精度只能保證前8 bit，所以，本文使用泰勒展開公式對精度進(jìn)行補(bǔ)償。最后，本文算法既保證了計算速度，又保證了精度。

本文在基于ARM架構(gòu)的鯤鵬920處理器上測試誤差結(jié)果,2個對比對象是本文算法和ARM V8架構(gòu)提供的求平方根函數(shù)的sqrt指令，結(jié)果如表3所示。

Table 3 Calculation accuracy

通過表3的測試誤差對比情況可以看出，其單精度浮點(diǎn)數(shù)計算結(jié)果的平均誤差為0.11E-08，而最大誤差為1.15E-07。雙精度浮點(diǎn)數(shù)計算結(jié)果的最大誤差為0.20E-15，且雙精度浮點(diǎn)數(shù)和雙精度復(fù)數(shù)的平均誤差為0.00E+00，最大誤差不超過4.40E-16。

7 結(jié)束語

本文介紹并總結(jié)了基于SIMD的平方根函數(shù)高性能實(shí)現(xiàn)和優(yōu)化原理與性能，并與已有的libm算法庫和處理器提供的sqrt指令進(jìn)行了比較，其總體性能是有優(yōu)勢的。針對平方根函數(shù)的實(shí)現(xiàn)主要是運(yùn)用了SIMD NEON語言，在ARM架構(gòu)上進(jìn)行并行化，之所以使用這項技術(shù)一方面是由于ARM架構(gòu)在市場上的需求逐漸提高，構(gòu)建ARM架構(gòu)完整的軟件生態(tài)十分重要，另一方面是由于目前SIMD技術(shù)的重要地位。所以，基于ARM架構(gòu)的函數(shù)開發(fā)是很有價值的。而今后的研究工作并不會只局限于本文中的平方根函數(shù)，接下來會繼續(xù)研究反三角函數(shù)、三角函數(shù)等其他函數(shù)的開發(fā)。