基于Piecewise算法的反正切運(yùn)算器的設(shè)計(jì)*

龍科蒞，汪 東，陳 虎，李旭軍

(1.湘潭大學(xué)物理與光電工程學(xué)院，湖南 湘潭 411105；2.湖南轂梁微電子有限公司，湖南 長(zhǎng)沙 410000)

1 引言

在科學(xué)計(jì)算、圖像處理和數(shù)字信號(hào)處理等應(yīng)用領(lǐng)域中，反正切運(yùn)算都是常用的基礎(chǔ)操作。由于軟件方法的運(yùn)算速度通常較慢，不能滿足應(yīng)用的實(shí)時(shí)性要求，需要設(shè)計(jì)專用于反正切運(yùn)算的硬件電路，即反正切運(yùn)算器。

在設(shè)計(jì)反正切運(yùn)算器的過(guò)程中，如何在運(yùn)算器的精度、速度和硬件成本三者之間進(jìn)行權(quán)衡是一個(gè)關(guān)鍵問(wèn)題。為了解決這一問(wèn)題，在以往的工作中提出了3類方法，即簡(jiǎn)單查表法[1,2]、數(shù)字迭代法[3 - 6]和Piecewise算法[7 - 13]。

簡(jiǎn)單查表法是指根據(jù)輸入從預(yù)先存儲(chǔ)的數(shù)據(jù)表中查找出對(duì)應(yīng)的運(yùn)算結(jié)果輸出。這種方法所需的表空間會(huì)隨所要求的運(yùn)算結(jié)果精度的提高而呈指數(shù)增長(zhǎng)。為了壓縮表空間，一類通過(guò)加法鏈接多個(gè)查找表從而獲得運(yùn)算結(jié)果的多部表法[1,2]應(yīng)運(yùn)而生，但這種方法若用于精度較高的運(yùn)算，所需的表空間仍然較大。

CORDIC(COordinate Rotation DIgital Computer)算法是典型的數(shù)字迭代法，其硬件耗費(fèi)較小,每次迭代只能獲得結(jié)果的一個(gè)數(shù)位，因此延時(shí)較高[3,4]。為了降低運(yùn)算延時(shí)，以增加運(yùn)算電路的復(fù)雜性為代價(jià)，Kebbati等人[5]開發(fā)了該算法的并行性；Lakshmi等人[6]則提高了CORDIC算法的運(yùn)算基數(shù)，使得一次迭代能獲得結(jié)果的更多數(shù)位。

Piecewise算法將運(yùn)算區(qū)域劃分成有限數(shù)量的子區(qū)間，在每一個(gè)子區(qū)間內(nèi)對(duì)原函數(shù)用多項(xiàng)式逼近以獲得誤差允許范圍內(nèi)的運(yùn)算結(jié)果。所有子區(qū)間內(nèi)用于逼近的多項(xiàng)式形式相同，僅根據(jù)子區(qū)間的不同來(lái)選取不同系數(shù)[7 - 9]。因此，子區(qū)間的數(shù)量直接決定了所需存儲(chǔ)的系數(shù)數(shù)量。另外，根據(jù)所劃分的子區(qū)間大小是否相同，又分為均勻[10,11]與非均勻[12,13]2種子區(qū)間劃分方式。非均勻劃分的方式更加靈活，往往能極大地減少子區(qū)間數(shù)量，從而壓縮系數(shù)表空間。Piecewise算法適合采用流水線的實(shí)現(xiàn)方式，且可通過(guò)調(diào)整分段子區(qū)間數(shù)量等方式來(lái)實(shí)現(xiàn)硬件成本與延時(shí)的平衡，具備較大的靈活性。

本文設(shè)計(jì)的反正切運(yùn)算器輸出結(jié)果與準(zhǔn)確結(jié)果之間的允許誤差為1 ulp (unit in the last place，最后位置單位)，足以滿足諸多實(shí)際應(yīng)用的要求。然而，在這種精度條件下，運(yùn)算器的輸出可能會(huì)隨輸入的遞增在可允許的誤差范圍內(nèi)變化，而不是像原函數(shù)那樣隨輸入的增加而逐步遞增，因此會(huì)引發(fā)單調(diào)性違例。在圖像處理等應(yīng)用中，運(yùn)算結(jié)果保持與原函數(shù)一致的單調(diào)性是十分重要的，比如紋理映射對(duì)于反正切運(yùn)算的輸出是否具備與原函數(shù)相一致的單調(diào)性十分敏感[14]。為了消除Piecewise算法引發(fā)的單調(diào)性違例，文獻(xiàn)[1]提出了一種通過(guò)調(diào)節(jié)斜率(Appropriate Slopes)和表初始值TIV(Table of Initial Values)來(lái)改善輸出的單調(diào)性方法；文獻(xiàn)[14]則通過(guò)測(cè)試所有輸出獲得引發(fā)單調(diào)違例的輸入范圍，并以此對(duì)系數(shù)表進(jìn)行經(jīng)驗(yàn)式的調(diào)整，以改善運(yùn)算的單調(diào)性。

Table 1 Subintervals obtained by two-level hierarchical segmentation表1 二級(jí)分層分段獲得的子區(qū)間

Table 2 Number of subintervals obstained by different segmentation methods表2 不同分段方法產(chǎn)生的子區(qū)間數(shù)量

Table 3 Signal width obtained by static analysis表3 靜態(tài)分析獲得的信號(hào)位寬

Table 4 Signal width obtained by dynamic analysis表4 動(dòng)態(tài)分析獲得的信號(hào)位寬

Table 5 Comparison results based on programmable logic表5 僅基于可編程邏輯的比較結(jié)果

本文基于Piecewise算法，設(shè)計(jì)了一種輸入輸出都符合IEEE-754單精度標(biāo)準(zhǔn)的浮點(diǎn)反正切運(yùn)算器，同時(shí)，使用二級(jí)分層分段方法進(jìn)行非均勻分段，以壓縮系數(shù)表空間。另外，通過(guò)控制單調(diào)性違例區(qū)域的尺度使其小于該區(qū)域內(nèi)輸入數(shù)的1 ulp大小，從而消除Piecewise逼近法帶來(lái)的單調(diào)性違例。此外，通過(guò)誤差分配以及對(duì)誤差與位寬之間關(guān)系的分析完成了位寬設(shè)計(jì)，并通過(guò)有限范圍內(nèi)的動(dòng)態(tài)分析進(jìn)行了位寬優(yōu)化。

2 算法與架構(gòu)

本文旨在完成如式(1)所示的反正切運(yùn)算。運(yùn)算的輸入輸出信號(hào)皆為IEEE-754標(biāo)準(zhǔn)的單精度浮點(diǎn)數(shù)。

(1)

2.1 多項(xiàng)式逼近

Piecewise算法往往在每個(gè)分段子區(qū)間內(nèi)采用一階、二階或三階多項(xiàng)式對(duì)目標(biāo)函數(shù)進(jìn)行可允許誤差范圍內(nèi)的逼近。由于“龍格效應(yīng)”，即逼近誤差不隨多項(xiàng)式的階數(shù)增加而一致減小，更高階的多項(xiàng)式很少被采用。根據(jù)文獻(xiàn)[8]可知，對(duì)于單精度浮點(diǎn)運(yùn)算，在硬件成本與延時(shí)之間進(jìn)行權(quán)衡，相較于需要一個(gè)較大系數(shù)表的Piecewise一階多項(xiàng)式逼近或者需要更多乘法和加法運(yùn)算的Piecewise三階多項(xiàng)式逼近，Piecewise二階多項(xiàng)式逼近更具優(yōu)勢(shì)。因此，本文在每個(gè)分段子區(qū)間內(nèi)使用二階多項(xiàng)式進(jìn)行運(yùn)算，以逼近該子區(qū)間內(nèi)的目標(biāo)函數(shù)。

假設(shè)f(x)=arctan(x)/(2π),x∈[a,b)。[a,b)被劃分為n個(gè)連續(xù)的子區(qū)間[ai,ai+1)，其中i=0,1,2,…,n-1。在每個(gè)子區(qū)間內(nèi)采用的二階多項(xiàng)式的統(tǒng)一形式如式(2)所示：

R=C0+D×(C1+D×C2)

(2)

其中,C0,C1和C2為多項(xiàng)式的系數(shù)，D是多項(xiàng)式的變量，R是多項(xiàng)式的運(yùn)算結(jié)果。式(2)所示的二階多項(xiàng)式形式相較于典型的二階多項(xiàng)式形式減少了一個(gè)乘法運(yùn)算。

假設(shè)xi是處于子區(qū)間[ai,ai+1)內(nèi)的一個(gè)給定值，對(duì)函數(shù)f(x)進(jìn)行二階泰勒級(jí)數(shù)展開，可以獲得式(2)中的系數(shù)和變量，如式(3)所示：

(3)

2.2 架構(gòu)

圖1展示了反正切運(yùn)算器的硬件架構(gòu)。圖1中壓縮單元根據(jù)反正切函數(shù)圖像的中心對(duì)稱性將輸入映射到[0,1]。地址單元可根據(jù)信號(hào)x的階碼以及尾數(shù)的高m位獲取xi,并通過(guò)譯碼獲得地址信號(hào)Addr。查找表單元可根據(jù)信號(hào)Addr分別讀取3個(gè)系數(shù)。減法單元完成D=x-xi的操作。之后2次乘加運(yùn)算完成如式(2)所示的多項(xiàng)式運(yùn)算。最后，規(guī)格化單元進(jìn)行就近舍入和規(guī)格化操作，以輸出符合IEEE-754標(biāo)準(zhǔn)的單精度浮點(diǎn)形式的結(jié)果。

Figure 1 Architecture of the arctangent operator圖1 反正切運(yùn)算器的硬件架構(gòu)

3 二級(jí)分層分段

3.1 逼近誤差分析

為了使反正切運(yùn)算器的輸出結(jié)果與標(biāo)準(zhǔn)結(jié)果之間的誤差不超過(guò)1 ulp，又因?yàn)楸疚脑O(shè)計(jì)對(duì)輸出結(jié)果使用就近舍入法最大可引入0.5 ulp的舍入誤差，因此逼近誤差必須小于0.5 ulp。根據(jù)泰勒定理，可以得到二階多項(xiàng)式逼近的誤差余項(xiàng)ε2，如式(4)所示：

(4)

其中，ξ是處于數(shù)軸上x到xi之間的一個(gè)數(shù)。

因此，逼近誤差εa可以由式(5)給出：

(5)

(6)

為了方便譯碼并減少子區(qū)間數(shù)量，本文使用二級(jí)分層分段方法，即在[2j,2j+1)內(nèi)，其中j為負(fù)整數(shù)，進(jìn)行均勻子區(qū)間劃分，并使xi=(ai+1-ai)/2，且子區(qū)間大小隨j的減小而逐步減小，以滿足精度要求。根據(jù)式(6)獲取每一個(gè)[2j,2j+1)內(nèi)滿足εa<0.5 ulp要求的子區(qū)間的長(zhǎng)度，由此可得如表1所示的子區(qū)間劃分。

特別地，對(duì)于子區(qū)間[0,2-12)內(nèi)設(shè)置xi=0，此時(shí)根據(jù)式(2)和式(3)可得R=x×(1/(2π))；對(duì)于子區(qū)間[2-12,2-11)，本文設(shè)置xi=2-12。此外，對(duì)于x=1，設(shè)置C2=C1=0,C0=0.125，此時(shí)輸出為0.125。

若使用均勻分段的方法，要達(dá)到相同精度，所需要的分段子區(qū)間的個(gè)數(shù)為1/2-12=212。表2展示了使用本文所述的二級(jí)分層分段進(jìn)行非均勻分段和使用均勻分段方法所分別產(chǎn)生的子區(qū)間個(gè)數(shù)。

3.2 單調(diào)性分析

圖2展示了使用如表1所示的二級(jí)分層分段子區(qū)間的Piecewise二階逼近算法在2個(gè)相鄰子區(qū)間[ai,ai+1)和[ai+1,ai+2)的交界處的單調(diào)性違例，且單調(diào)性違例區(qū)域大小為α=α1+α2。在ai+1的某個(gè)鄰域內(nèi)，任意相鄰2個(gè)單精度浮點(diǎn)形式的數(shù)的距離為β=ULP(ai+1)，其中ai+1為單精度浮點(diǎn)形式。由圖2可知，當(dāng)β≥α?xí)r，單調(diào)性違例區(qū)域?qū)⒉粫?huì)對(duì)輸出的單調(diào)性產(chǎn)生影響。

Figure 2 Monotonic violation at the junction of adjacent subintervals圖2 2個(gè)子區(qū)間交界處的單調(diào)性違例

為了獲得α2，假設(shè)：

(7)

如圖2所示，當(dāng)f(x1)=f(x2)時(shí)，根據(jù)式(2)和式(3)并忽略α2的高次項(xiàng)可得式(8)：

(8)

其中

(9)

(10)

(11)

4 位寬設(shè)計(jì)與優(yōu)化

4.1 靜態(tài)分析

除舍入誤差和逼近誤差外，設(shè)計(jì)還需要考慮由于信號(hào)位寬有限而引入的量化誤差。根據(jù)表1以及式(5)和式(6)可以得到逼近誤差εa≤2-4ulp；又由于可允許的最大誤差不能超過(guò)1 ulp，輸出的最大舍入誤差為0.5 ulp,因此必須使得量化誤差εq≤(2-2+2-3+2-4) ulp。在如圖1展示的架構(gòu)中，量化誤差來(lái)自于信號(hào)C0,C1,C2以及M0,M1,M2，假設(shè)其值分別為εC0,εC1,εC2以及εM0,εM1,εM2。根據(jù)圖1、表1、式(2)和式(3)，可以得到式(12)：

(12)

根據(jù)式(3)可以得到22C2≥f(x)，又由于f(x)≤x/(2π)，則22ULP(C2)≥ULP(f(x))，其中f(x)為IEEE-754標(biāo)準(zhǔn)的單精度浮點(diǎn)數(shù)；更進(jìn)一步地，根據(jù)式(2)可以得到εC2≤22ULP(C2)×D2；又根據(jù)表1可知,本文所述分段方式下D≤2-8-1，因此εC2≤ULP(C2)×2-16。

在位寬設(shè)計(jì)過(guò)程中，為了減少硬件成本且避免不必要的精度損失，本文分配了較大的量化誤差給遠(yuǎn)離輸出端的信號(hào)，如圖1中所示的信號(hào)C2和M2；同時(shí)，分配較小的量化誤差給遠(yuǎn)離輸入端的信號(hào)，如信號(hào)C0和M0；此外，分配更大的量化誤差給乘法運(yùn)算輸入，而分配更小的量化誤差給加法運(yùn)算輸入。根據(jù)以上量化誤差分配原則，本文分配的量化誤差如式(13)所示：

(13)

由此，根據(jù)式(12)可以得到如表3所示的信號(hào)位寬。根據(jù)式(12)可知εC2≤ULP(C2)×2-16，又由于εC2=2-3ulp，則ULP(C2)≥2-13。因此，信號(hào)C2的尾數(shù)位數(shù)為24-13=11 bit，其中算式中的24是指IEEE-754標(biāo)準(zhǔn)的單精度浮點(diǎn)形式的輸出具備24位尾數(shù)。因此，信號(hào)C2的位寬為1+8+11=20 bit，其中包含1 bit符號(hào)位，8 bit階碼位。

4.2 動(dòng)態(tài)分析

為了優(yōu)化信號(hào)位寬以壓縮硬件成本，本文對(duì)[2-12,1]內(nèi)的所有單精度浮點(diǎn)形式的輸入以及對(duì)應(yīng)的輸出進(jìn)行動(dòng)態(tài)分析，其主要步驟如下所示：

步驟1減少指定信號(hào)位寬。

步驟2測(cè)試輸出結(jié)果與黃金模型運(yùn)算結(jié)果之間的誤差是否不大于1 ulp，以及輸出是否隨輸入的增加而逐步遞增。如果測(cè)試通過(guò)則位寬減少成功；否則將信號(hào)位寬恢復(fù)到步驟1之前的狀態(tài)。

步驟3重復(fù)步驟1和步驟2，直到完成所有指定信號(hào)的位寬優(yōu)化。

5 仿真結(jié)果

由于反正切函數(shù)圖像具備中心對(duì)稱性，且在[0,2-12)內(nèi)Output=Input×(1/(2π))。因此，對(duì)于[2-12,1]內(nèi)的所有輸入(Input)進(jìn)行的測(cè)試被認(rèn)為具備覆蓋性，測(cè)試結(jié)果如圖3所示。誤差Error=|Output-Std|/ULP(Std)，其中Std為C語(yǔ)言數(shù)學(xué)庫(kù)中對(duì)應(yīng)函數(shù)的運(yùn)算結(jié)果。由圖3可知，反正切運(yùn)算器的輸出與黃金模型運(yùn)算結(jié)果之間的誤差不超過(guò)1 ulp，達(dá)到了設(shè)計(jì)精度要求。

Figure 3 Error of the output of all inputs within [2-12,1]圖3 [2-12,1]內(nèi)所有輸入的輸出結(jié)果的誤差

類似地，本文對(duì)[2-13,1]內(nèi)的所有輸入進(jìn)行單調(diào)性驗(yàn)證。驗(yàn)證方法是將2個(gè)相鄰輸入(Input)的輸出結(jié)果(Output)相減得到差(Differences)，其中對(duì)應(yīng)于較大輸入的輸出作為被減數(shù)，獲得的結(jié)果如圖4所示。根據(jù)圖4可知，所有的差都不小于0，即輸出結(jié)果隨著輸入的增大而逐步增大，具備與原函數(shù)一致的單調(diào)性，達(dá)到了設(shè)計(jì)要求。

Figure 4 Monotonicity test results of the output corresponding to all inputs within [2-13,1]圖4 [2-13,1]內(nèi)所有輸入對(duì)應(yīng)的輸出的單調(diào)性測(cè)試結(jié)果

進(jìn)一步地，本文展示了所設(shè)計(jì)的反正切運(yùn)算器與芯片TMS320F2837xD的指令A(yù)TANPUF32在輸入處于(0.031017+0.3×10-7,0.031017+1.6×10-7)內(nèi)的輸出隨輸入變化的圖像，分別如圖5和圖6所示。從圖5可以知道，TMS320F2837xD的指令A(yù)TANPUF32的輸出并不隨輸入的增加而逐步地遞增，即存在單調(diào)性違例。從圖6可知，在同樣的輸入范圍內(nèi)，本文所設(shè)計(jì)的反正切運(yùn)算器的輸出隨輸入的增加而逐步遞增，且輸出隨輸入變化的曲線相對(duì)于圖5更加平滑。

Figure 5 Output of the instruction ATANPUF32 in the TMS320F2837xD changing with the input圖5 TMS320F2837xD的指令A(yù)TANPUF32的輸出隨輸入變化的圖像

Figure 6 Output of the designed arctangent unit changing with the input圖6 本文設(shè)計(jì)的反正切運(yùn)算器的輸出隨輸入變化的圖像

為了比較硬件成本，本文使用XILINX Virtex-4 XC4VLX100-12 FPGA進(jìn)行了硬件仿真。此外，為了獲得無(wú)偏見的比較結(jié)果，該硬件仿真僅僅基于可編程的邏輯單元。本文還與文獻(xiàn)[15]的設(shè)計(jì)進(jìn)行了比較，其使用了與本文類似的Piecewise二階逼近算法和二次乘加運(yùn)算結(jié)構(gòu)來(lái)完成類似的非線性函數(shù)運(yùn)算。如表5所示，本文設(shè)計(jì)的反正切運(yùn)算器在面積和延時(shí)2個(gè)方面都優(yōu)于文獻(xiàn)[15]的類似設(shè)計(jì)。

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

本文基于Piecewise算法設(shè)計(jì)了一個(gè)輸入輸出都符合IEEE-754單精度浮點(diǎn)標(biāo)準(zhǔn)的反正切運(yùn)算器。逼近誤差分析與二級(jí)分層分段方法相結(jié)合完成了運(yùn)算區(qū)域的非均勻分段，達(dá)到了壓縮系數(shù)表的目的。通過(guò)分析單調(diào)性違例區(qū)域的尺寸與對(duì)應(yīng)子區(qū)間大小之間的關(guān)系，使得單調(diào)性違例區(qū)域的尺寸可控，從算法上確保了運(yùn)算結(jié)果的單調(diào)性。此外，通過(guò)分析量化誤差完成了關(guān)鍵路徑上信號(hào)位寬的設(shè)計(jì)，并進(jìn)一步通過(guò)有限范圍內(nèi)動(dòng)態(tài)誤差分析完成了位寬優(yōu)化，達(dá)到了減少硬件成本的目的。仿真驗(yàn)證結(jié)果表明，本文設(shè)計(jì)的反正切運(yùn)算器的硬件開銷相較于同類設(shè)計(jì)更小，輸出結(jié)果與黃金模型運(yùn)算結(jié)果之間的誤差小于1 ulp，且具備與原函數(shù)一致的單調(diào)性。