仲雅莉 ，吳俊輝 ，劉炫 ，高萍 ,3，段曉輝 ，4

(1.國家超級(jí)計(jì)算無錫中心，江蘇 無錫 214072；2;江南大學(xué)，江蘇 無錫 214122；3.山東大學(xué)，山東 濟(jì)南 250100；4.清華大學(xué)，北京 100084)

0 引言

坐標(biāo)旋轉(zhuǎn)計(jì)算機(jī)(Coordinated Rotation Digital Computer，CORDIC）算法只有移位和加減運(yùn)算，便于在FPGA 等硬件平臺(tái)上實(shí)現(xiàn)復(fù)雜的三角函數(shù)、雙曲函數(shù)、指數(shù)函數(shù)和復(fù)數(shù)求模等計(jì)算，廣泛應(yīng)用于數(shù)字鑒相、數(shù)字上下變頻、波形產(chǎn)生、快速傅里葉變換等方面[1-4]。

針對(duì)CORDIC 算法中由于迭代次數(shù)多、輸出延時(shí)大和精度較低等問題，國內(nèi)外很多學(xué)者進(jìn)行了研究和改進(jìn)。文獻(xiàn)[5]提出了基于自適應(yīng)旋轉(zhuǎn)角度的CORDIC 算法，該算法雖然減少了迭代次數(shù)，但每次迭代都需要額外判斷，增加了實(shí)現(xiàn)難度。文獻(xiàn)[6]提出了將查找表和傳統(tǒng)CORDIC 算法相融合，通過查找表將角度值細(xì)化，通過數(shù)學(xué)量化分析，根據(jù)細(xì)化后的較小角度補(bǔ)碼，直接按位值進(jìn)行角度單向旋轉(zhuǎn)，該算法隨著對(duì)精度要求的提高，查找表存儲(chǔ)空間大大增加，硬件資源消耗巨大。文獻(xiàn)[7]提出了一種基于最佳一致逼近方法的幅度與相位補(bǔ)償算法，第一步是利用傳統(tǒng)的CORDIC 算法迭代數(shù)次后得到向量信息，第二步是采用最佳逼近法進(jìn)行多項(xiàng)式補(bǔ)償，該算法雖然提高了精度，但逼近算法需要存儲(chǔ)多項(xiàng)式系數(shù)，硬件資源消耗較大。文獻(xiàn)[8]提出來一種基于CORDIC 算法的反正切函數(shù)計(jì)算的改進(jìn)算法，該算法對(duì)累加器中因截尾而產(chǎn)生的誤差做了算法改進(jìn)，僅增加了運(yùn)算速度，精度并未明顯提高。文獻(xiàn)[9]提出了一種基于查找表的改進(jìn)的CORDIC 算法，該方法通過縮減有效數(shù)據(jù)位寬、合并迭代等手段節(jié)省了剩余角度Z 的計(jì)算量，該算法把迭代次數(shù)進(jìn)行拆分，用了14 次迭代和查找表的結(jié)合，硬件資源消耗較多，性能提升卻不足。

在分子動(dòng)力學(xué)模擬中[10-11]，成鍵力的計(jì)算引入了三角函數(shù)。常見的一種情況是從三維坐標(biāo)計(jì)算角度，通常需要根據(jù)sinθ 和cosθ 的值求解反正切，該應(yīng)用場景下對(duì)反正切的精度要求較高，傳統(tǒng)的CORDIC 算法很難滿足[12-13]。

本文在文獻(xiàn)[7]的基礎(chǔ)上，簡化逼近算法，提出了一種改進(jìn)算法，測得不采用函數(shù)擬合也可得到高精度的剩余相位跟蹤，從而有效提高了CORDIC 求反正切算法的精度，降低了延時(shí)，減少了電路面積。

1 CORDIC 算法基本原理

在X-Y 坐標(biāo)平面上將向量(x0，y0)逆向旋轉(zhuǎn)角度θ 得到向量(x1，y1)，如圖1 所示。兩向量間坐標(biāo)變換關(guān)系為：

圖1 向量旋轉(zhuǎn)示意圖

兩邊同時(shí)除以cosθ，可得到偽旋轉(zhuǎn)方程，即：

除以cosθ 后，旋轉(zhuǎn)角度還是正確的，但向量的模長發(fā)生了變化。

CORDIC 求反正切的本質(zhì)就是從坐標(biāo)(x0，y0)開始旋轉(zhuǎn)，每次旋轉(zhuǎn)固定角度θi，旋轉(zhuǎn)方向?yàn)閐i=sign(y(i))，旋轉(zhuǎn)的目標(biāo)是縱坐標(biāo)yi趨近于0，進(jìn)行N 步迭代，則有反正切角。為了能夠簡化運(yùn)算部件，CORDIC算法常采用每一小步的旋轉(zhuǎn)角度θi滿足tanθi=2-i，從而可以通過簡單移位來完成乘法運(yùn)算，因此原始的算法可以簡化為迭代移位相加算法，迭代方程為：

如果需要保持向量的模長，則可以通過加入旋轉(zhuǎn)補(bǔ)償因子K 實(shí)現(xiàn)。當(dāng)N 確定時(shí)，經(jīng)歷N 次迭代后的K 是一個(gè)常數(shù)：

該補(bǔ)償因子常在求三角函數(shù)sin 和cos 時(shí)使用。

2 基于CORDIC 求反正切的算法改進(jìn)

對(duì)CORDIC 算法進(jìn)行分析可以看出，θi=arctan2-i。而進(jìn)行i 次迭代后，殘余的角度的絕對(duì)值即是求反正切值時(shí)的誤差。定義εi為第i 次旋轉(zhuǎn)后的誤差，則：

所以CORDIC 求反正切算法在殘余角度較小時(shí)的誤差以每次1 bit 的速度收斂。對(duì)于定點(diǎn)數(shù)而言，即每次迭代使得精度收斂1 bit。以在進(jìn)行均勻取樣仿真了不同迭代次數(shù)下的CORDIC 最大誤差如表1 所示，對(duì)于相位精度要求極高的場景，如果要達(dá)到10-9數(shù)量級(jí)，需要迭代28 次，其帶來的硬件資源消耗和時(shí)延都是不可承受的。

表1 不同迭代次數(shù)下CORDIC 算法的最大誤差

2.1 使用sinε 逼近ε

當(dāng)CORDIC 經(jīng)過第N 次旋轉(zhuǎn)之后，設(shè)其估計(jì)的角度剩余誤差為ε，顯然：

即N 次迭代后的殘余誤差小于2-(N-1)。根據(jù)微積分中正弦函數(shù)的極限，有：

所以使用sinε 來逼近ε 是一種有潛力降低誤差的手段。

由于sinε 與ε 在ε 比較小時(shí)非常接近，一種可以考慮的方案是直接采用sinε 來代替ε。該方案的誤差可以通過泰勒展開進(jìn)行估計(jì)，ε-sinε 在ε=0 處展開可得：

也就是說，sinε 和ε 的差值以三次方的速度縮小，比直接CORDIC的收斂速度快得多。當(dāng)ε接近2-9≈1.95×10-3時(shí)，從式(9)可以估計(jì)出ε-sinε≈1.24×10-9，意味著用CORDIC算法迭代10次后，反正切的誤差約為1.95×10-3，而使用sinε 對(duì)反正切結(jié)果進(jìn)行校正后，誤差降為約1.24×10-9。

實(shí)驗(yàn)的結(jié)果可以與上述估計(jì)相吻合：令err=ε-sinε，畫出曲線如圖2所示，可以看出ε越小，ε與sinε的差值越小。當(dāng)ε＜2-9時(shí)，ε與sinε的差值可達(dá)到10-9數(shù)量級(jí)，這也意味著當(dāng)CORDIC迭代10次后，殘余的ε足夠小，可以用sinε的值代替ε的值來修正CORDIC 的結(jié)果。

圖2 區(qū)間[2-14，2-3]上的ε-sinε(對(duì)數(shù)坐標(biāo)軸)

根據(jù)上述對(duì)誤差的論證，在一定迭代步數(shù)后采用sinε 代替ε 是安全的，如果輸入的初始向量已經(jīng)在單位圓上，那么最終向量的縱坐標(biāo)即為sinε。所以本文對(duì)應(yīng)修改了CORDIC 求反正切的算法，如算法1 所示。

算法1輸入在x 軸正半軸單位圓上的改進(jìn)CORDIC算法

在迭代步數(shù)N 的選取上，可以根據(jù)式(7)和式(9)，按照所需的誤差t 來對(duì)N 進(jìn)行選?。?/p>

而對(duì)于普通的CORDIC 算法，第N 步的誤差約為2-(N-1)，則需要：

從式(11)和式(12)可以看出，在需求精度較高的情況下，本文改進(jìn)的CORDIC 算法只需要大約三分之一的迭代步就能實(shí)現(xiàn)與原有的CORDIC 相似的精度。

2.2 sinε 的計(jì)算

根據(jù)CORDIC 的原理可知，CORDIC 的旋轉(zhuǎn)過程中向量模長會(huì)發(fā)生變化，需要對(duì)模長進(jìn)行修正。設(shè)迭代初始向量為(x0，y0)，經(jīng)過N 步旋轉(zhuǎn)后，根據(jù)式(3)得到(xN，yN)，zN，式(4)可求得旋轉(zhuǎn)修正因子K 使得：

當(dāng)(x0，y0)在單位圓上時(shí)，如圖3 所示，l=1，yN即為sinε，這就可以解決輸入的坐標(biāo)在單位圓上的問題。

圖3 (x0，y0)在單位圓上時(shí)的余角

若輸入的坐標(biāo)(x0，y0)不在單位圓上，則需要根據(jù)式(14)進(jìn)行模值歸一化：

從而使得：

如式(14)所示，模值歸一化引入了平方根倒數(shù)，直接求解資源消耗大。可以用牛頓迭代法求得，迭代公式如式(17)所示，其中p 是式(14)中初始坐標(biāo)(x0，y0)的模值平方，q 是對(duì)其平方根倒數(shù)的估計(jì)。

在定點(diǎn)數(shù)實(shí)現(xiàn)牛頓迭代法時(shí)，由于沒有尾碼對(duì)齊的性質(zhì)[14]，因此很難通過簡單的變換就獲得一個(gè)良好的初值[15]。本文考慮將p 映射到[1，4)之間，如式(18)所示：

其中22k和2k均可通過移位實(shí)現(xiàn)。當(dāng)p'落在[1，4)之間時(shí)，可以分段通過線性變換獲得一個(gè)較為良好的初值：

式(19)中初值與結(jié)果的誤差基本可以控制在2×10-2以內(nèi)，如圖4 所示。這樣只需兩步牛頓迭代就可以獲得較為準(zhǔn)確的。由于CORDIC 旋轉(zhuǎn)與模值無依賴性，模值歸一化流程僅在最后一步將s 與Kyn相乘時(shí)會(huì)導(dǎo)致額外的一次乘法的時(shí)延。

圖4 初值q0與之差

2.3 整體架構(gòu)

改進(jìn)的CORDIC 求反正切算法整體架構(gòu)如圖5 所示。其中，角度預(yù)處理是指把坐標(biāo)轉(zhuǎn)換到第一或第四象限，預(yù)處理后的坐標(biāo)先按照傳統(tǒng)CORDIC 進(jìn)行N 級(jí)迭代，再根據(jù)式(16)進(jìn)行縱坐標(biāo)修正和誤差補(bǔ)償，最后結(jié)果映射還原到原始象限的角度。算法2 展示了該架構(gòu)對(duì)應(yīng)的偽代碼。

圖5 改進(jìn)的CORDIC 算法的整體流程

算法2整體的CORDIC 算法架構(gòu)

2.4 改進(jìn)算法在分子動(dòng)力學(xué)模擬中的應(yīng)用

在分子動(dòng)力學(xué)模擬中，常見的一種情況是從三維坐標(biāo)計(jì)算角度。以分子動(dòng)力學(xué)中的二面角? 為例，很多實(shí)現(xiàn)就選擇了使用sin? 和cos? 作為二維反正切函數(shù)的輸入，從而獲得值域?yàn)閇-π，π)的結(jié)果。此時(shí)由于坐標(biāo)顯然在單位圓上，模長歸一化的過程就可以省略。

圖6 展示了求解平面ijk 和jkl 間二面角? 所需的向量，r 是兩個(gè)原子間的向量，n 是平面的法向量。其中：

圖6 二面角求解時(shí)的向量

n'是nijk和njkl所確定的平面上與nijk垂直的向量：

其二維視圖如圖7 所示。這樣，? 的計(jì)算方式為：

圖7 圖6 中幾個(gè)向量的二維視圖

也就是說此處的模長歸一化已經(jīng)隱含在了二面角計(jì)算的流程中，從而可以省略反正切函數(shù)的模長歸一化過程。

類似地，在CHARMM 鍵角的計(jì)算中：

這里V是鍵角的勢能，K和θ0是經(jīng)驗(yàn)參數(shù)，θ是兩個(gè)鍵形成鍵角的角度，對(duì)于以i 原子為中心，j、k 構(gòu)成的鍵角，一般以反余弦函數(shù)求解：

計(jì)算原子受力F 時(shí)需要根據(jù)鏈?zhǔn)椒▌t求解V 對(duì)三個(gè)原子坐標(biāo)x 的梯度，以j 原子為例：

也就是說，cosθ 本來就需要計(jì)算，通過對(duì)計(jì)算流程進(jìn)行修改后即可避免反正切函數(shù)進(jìn)行模長歸一化的過程。

綜上所述，雖然在必要時(shí)可以利用牛頓迭代法實(shí)現(xiàn)模長歸一化的過程，但是在分子動(dòng)力學(xué)等求解具體計(jì)算幾何問題的情境下，基本不需要額外的模長歸一化實(shí)現(xiàn)。此時(shí)，改進(jìn)的CORIDC 算法可以更好地發(fā)揮其優(yōu)勢。

3 算法仿真和分析

將傳統(tǒng)CORDIC 求反正切算法和改進(jìn)的CORDIC 求反正切算法分別用MATLAB 建模仿真，并與理論值進(jìn)行了誤差分析。輸入數(shù)據(jù)的橫坐標(biāo)和縱坐標(biāo)分別以0.001為間隔遍歷[-1，1]，統(tǒng)計(jì)相同迭代次數(shù)下所有數(shù)據(jù)的最大絕對(duì)誤差，性能曲線如圖8 所示?？梢钥闯龈倪M(jìn)后的CORDIC 算法性能提升明顯，在迭代10 次的情況下，精度從10-3數(shù)量級(jí)提升到10-9數(shù)量級(jí)。

圖8 改進(jìn)CORDIC 算法與傳統(tǒng)CORDIC 算法的誤差比較

本研究采用Verilog HDL 語言進(jìn)行設(shè)計(jì)，并在Xilinx公司的Virtex UltraScale+HBM VCU128 型號(hào)FPGA 上完成了電路實(shí)現(xiàn)。整個(gè)模塊的設(shè)計(jì)采用了32 bit(3 bit 整數(shù)位+29 bit 小數(shù)位)的定點(diǎn)化輸入輸出。在保證絕對(duì)誤差小于5×10-9的前提下，比較了硬件資源消耗如表2 所示。

表2 改進(jìn)CORDIC 算法與傳統(tǒng)CORDIC 算法的資源消耗對(duì)比

對(duì)于需要模值歸一化的應(yīng)用場景，歸一化處理引入了DSP 資源的消耗，但查找表(Look Up Table，LUT)資源消耗降低了64.8%，觸發(fā)器(Flip-Flop，F(xiàn)F)資源消耗降低了35.3%，輸出時(shí)延降低了53.3%。對(duì)于無模值歸一化應(yīng)用場景，LUT 資源消耗降低了64.3%，F(xiàn)F 資源消耗降低了61.6%，輸出時(shí)延降低了60%。改進(jìn)的CORDIC 算法在硬件消耗和時(shí)延方面有明顯優(yōu)勢。

4 結(jié)論

針對(duì)傳統(tǒng)CORDIC 求反正切函數(shù)硬件資源消耗大，迭代次數(shù)多等問題，實(shí)現(xiàn)了一種改進(jìn)的CORDIC 算法計(jì)算反正切函數(shù)的解法。在保證與傳統(tǒng)CORDIC 算法精度數(shù)量級(jí)相同的情況下，減少了算法迭代次數(shù)，用旋轉(zhuǎn)后的縱坐標(biāo)補(bǔ)償了相位，從而顯著提高了算法運(yùn)算速度。采用32 bit 定點(diǎn)的硬件實(shí)現(xiàn)，在保證相同計(jì)算精度的前提下，改進(jìn)CORDIC 算法的硬件資源消耗和時(shí)延方面均有明顯降低。算法仿真和FPGA 實(shí)際測試結(jié)果表明，改進(jìn)后的CORDIC 算法適用于高速高精度場合。

該設(shè)計(jì)在包含模長歸一化實(shí)現(xiàn)時(shí)已經(jīng)比經(jīng)典的CORDIC 算法擁有更短的延時(shí)，在很多計(jì)算幾何的具體問題中，可以通過對(duì)計(jì)算流程的適當(dāng)重構(gòu)從而使用單位向量作為輸入，從而可以進(jìn)一步降低該設(shè)計(jì)的開銷，尤其適合分子動(dòng)力學(xué)等三維空間中需要使用計(jì)算幾何方法求解的問題。