趙永浩,賈海鵬,張?jiān)迫?張思佳

(1.大連海洋大學(xué)信息工程學(xué)院，遼寧 大連 116023;2.中國(guó)科學(xué)院計(jì)算技術(shù)研究所計(jì)算機(jī)體系結(jié)構(gòu)國(guó)家重點(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ì)算機(jī)圖形學(xué)和積分計(jì)算等科學(xué)領(lǐng)域。隨著處理器架構(gòu)的不斷發(fā)展和完善，為了更好地支持和完善計(jì)算機(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ù)是對(duì)ARM Cortex-A和Cortex-R系列處理器中SIMD的擴(kuò)展。NEON寄存器存儲(chǔ)的是同一數(shù)據(jù)類型的向量，相比其他指令其優(yōu)點(diǎn)是NEON指令可以對(duì)4個(gè)32 bit和2個(gè)64 bit的元素同時(shí)進(jìn)行計(jì)算，并且還支持int、float、double等多種數(shù)據(jù)類型。該技術(shù)通過(guò)同時(shí)計(jì)算多個(gè)元素，加快計(jì)算速度，來(lái)達(dá)到提升整體性能的目的。

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

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

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

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

2 相關(guān)工作

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

2.1 libm

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

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

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

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

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

2.2 ARMPL

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

2.3 ARM V8 架構(gòu)

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

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

2.4 向量化SIMD

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

2.5 平方根算法綜述

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

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

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

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

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

3 研究背景

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

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

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

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

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

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

(1)

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

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

(2)

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

3.2 泰勒級(jí)數(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é)果很小時(shí)，則有：

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

(4)

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

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

(5)

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

(6)

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

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

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

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

(7)

(8)

(9)

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

(10)

最終表達(dá)式為：

(11)

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

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

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。

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

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

(12)

具體包括如下步驟：

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

(13)

(14)

最終的計(jì)算公式為：

(15)

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

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

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

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

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

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

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

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

(16)

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

(17)

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

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

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

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

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

算法2x0計(jì)算過(guò)程

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。

在計(jì)算x1時(shí)使用的是泰勒公式展開(kāi)的方法，具體計(jì)算過(guò)程如下所示:

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

x1=1-x*y2

(18)

(19)

具體計(jì)算過(guò)程如算法3所示。

算法3x1計(jì)算過(guò)程

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；

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

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

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

(20)

(21)

通過(guò)2個(gè)公式可以分別求出開(kāi)根號(hào)之后的實(shí)部和虛部結(jié)果。具體計(jì)算過(guò)程如算法4所示。

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

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ù)的計(jì)算時(shí)，可以利用ARM NEON INTRINSIC指令對(duì)SIMD代碼進(jìn)行優(yōu)化，以此加快對(duì)平方根函數(shù)的計(jì)算。此外，由于ARM NEON INTRINSIC提供的SIMD乘加指令的速度很快，因此可以在對(duì)程序進(jìn)行精度補(bǔ)償時(shí)使用泰勒展開(kāi)公式，從而使用多個(gè)SIMD乘加指令加速運(yùn)算。

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

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

6 性能評(píng)估

6.1 測(cè)試環(huán)境

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

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

Table 1 Hardware environment of experiment

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

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

6.2 性能分析

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

Table 2 Performance comparison

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

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

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

Table 3 Calculation accuracy

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

7 結(jié)束語(yǔ)

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