基于SIMD的Square Root函數(shù)高性能實現(xiàn)與優(yōu)化*

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

(1.大連海洋大學信息工程學院，遼寧 大連 116023;2.中國科學院計算技術研究所計算機體系結構國家重點實驗室，北京100190)

1 引言

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

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

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

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

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

2 相關工作

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

2.1 libm

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

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

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

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

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

2.2 ARMPL

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

2.3 ARM V8 架構

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

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

2.4 向量化SIMD

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

2.5 平方根算法綜述

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

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

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

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

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

3 研究背景

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

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

IEEE 754標準浮點數(shù)在內存中的存儲方式如圖1所示。

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

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

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

(1)

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

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

(2)

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

3.2 泰勒級數(shù)

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

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

(3)

當|x-x0|的結果很小時，則有：

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

(4)

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

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

(5)

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

(6)

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

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

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

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

(7)

(8)

(9)

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

(10)

最終表達式為：

(11)

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

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

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

Output:y。

1．a←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。

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

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

(12)

具體包括如下步驟：

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

(13)

(14)

最終的計算公式為：

(15)

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

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

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

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

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

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

通過以上運算，得出e的表達式為：

(1)當e% 2=0時：

(16)

(2)當e% 2=1時：

(17)

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

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

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

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

(4)平方根倒數(shù)算法得出的結果為float類型，需要將其轉換成double類型，可使用SIMD NEON提供的verinterpretq_f64_u64轉換指令完成。這條指令在轉換后沒有改變值，只改變數(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)在取出雙精度浮點數(shù)的尾數(shù)位的運算中，使用0x001fffffffffffff和64位無符號整數(shù)進行“與”運算,得到的結果會保留指數(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；

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

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

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

(20)

(21)

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

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

Input：z=a+bi。

Output：w=c+di。

1．a=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ù)算法的實現(xiàn)，在執(zhí)行平方根函數(shù)的計算時，可以利用ARM NEON INTRINSIC指令對SIMD代碼進行優(yōu)化，以此加快對平方根函數(shù)的計算。此外，由于ARM NEON INTRINSIC提供的SIMD乘加指令的速度很快，因此可以在對程序進行精度補償時使用泰勒展開公式，從而使用多個SIMD乘加指令加速運算。

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

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

6 性能評估

6.1 測試環(huán)境

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

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

Table 1 Hardware environment of experiment

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

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

6.2 性能分析

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

Table 2 Performance comparison

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

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

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

Table 3 Calculation accuracy

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

7 結束語

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