CORDIC

CORDIC（Coordinate rotation digital computer），也稱為Volder演演算法（Volder's algorithm），是一個可以計算三角函數，簡單且有效率的演演算法，可以在任意進制下運算，一般會每次計算一位數字。因此CORDIC屬於逐位計算（Digit-by-digit）方法中的一個例子。

CORDIC演演算法還有其他的名稱，像是圓形CORDIC (Jack E. Volder)、線性CORDIC、雙曲線CORDIC(John Stephen Walther)、及泛用雙曲線CORDIC（GH CORDIC, Yuanyong Luo et al.）。用類似的方式也可以計算雙曲函數、平方根、乘法、除法、指數及對數等。

CORDIC和一些名為「偽乘法」（pseudo-multiplication）、「偽除法」（pseudo-division）及factor combining的方法，常用在沒有乘法器的應用（像是簡單的微控制器以及FPGA），其中會用到的運算是加法、減法、位元移位及查找表。這些演算法可歸類在「移位和相加」（shift-and-add）演演算法中。在電腦科學中，若CPU沒有硬體的乘法器，常會用CORDIC來實現浮點數運算。

## 歷史

英國數學家亨利·布里格斯早在1624年時就已發現此演算法，後來Robert Flower也在1771年時發現，不過後來是因為低複雜度的有限狀態CPU，才針對CORDIC作較進一步的最佳化。

CORDIC是在1956年問世，是由康維爾空氣電子部門的Jack E. Volder發現，一開始是因為要取代B-58轟炸機導航電腦上面的類比式解角器（resolver），更換成更準確而實時的數位方案。因此，有時也會將CORDIC稱為數位解角器（digital resolver）。

Volder的研究，是因為1946年版CRC化學和物理手冊中的公式而得到靈感：

$$\begin{align} <pre> K_n R \sin(\theta \pm \varphi) &= R \sin(\theta) \pm 2^{-n} R \cos(\theta), \\ K_n R \cos(\theta \pm \varphi) &= R \cos(\theta) \mp 2^{-n} R \sin(\theta), \\ </pre> \end{align}$$

其中$\varphi$是使$\tan(\varphi) = 2^{-n}$成立的值，且$K_n := \sqrt{1 + 2^{-2n}}$。

他的研究最後產生了一個內部的技術報告，提到用CORDIC演演算法來求解正弦及餘弦函數，以及一個實現此功能的原型電腦。報告中也提到用修改版的CORDIC演演算法計算雙曲函數、座標旋轉、對數及指數的可能性。用CORDIC來進行乘法和除法運算的想法也是在此時形成。依照CORDIC演演算法的原則，Volder的同事Dan H. Daggett發展了在二進位以及二進碼十進數（BCD）之間轉換的演演算法。

## 應用

CORDIC用簡單的移位和相加運算來處理像是三角函數、雙曲函數、對數函數、實數及複數乘法、除法、方根計算、線性系統求解、特徵值估測、奇異值分解、QR分解等。因此，CORDIC可以用在許多的應用中，像是信號處理、影像處理、通訊系統、機器人學及三維電腦圖形等。

### 硬體

若沒有硬體乘法器的話，CORDIC一般會比其他演算法要快很多，若是用FPGA或ASIC的話，使用的邏輯閘也會少很多。

CORDIC是FPGA開發應用程式（像是Xilinx的Vivado）中的標準半導體IP核，而不是使用特殊函數的冪級數實現，其原因是CORDIC IP的通用性，CORDIC可以計算許多不同的函數，而為特定函數開發的乘法器只能計算特定的函數。

另一方面，若有硬體乘法器（例如DSP），查表法及冪級數會比CORDIC快很多。近年來，CORDIC演演算法常用在生醫應用中，尤其是用FPGA實現的應用。

使用泰勒級數的問題是此方法可以產生小的絕對誤差，但其中沒有良態的相對誤差。其他多項式近似法，例如Minimax近似演演算法，可以同時控制這二種的誤差。

### 軟體

在CPU只有整數運算的古老系統中，會將CORDIC放在其IEEE 754函式庫的一部份。現代的通用CPU已有浮點運算暫存器，也有加法、減法、乘法、除法、三角函數、平方根、一般對數、自然對數等，幾乎沒有用到CORDIC的場合。只有一些有特殊安全性或是時間要求的應用程式才會用到CORDIC。

## 運作模式

### 旋轉模式

CORDIC可以用來計算許多不同的函數。以下說明如何在旋轉模式（rotation mode）下的CORDIC計算角度的正弦函數和餘弦函數，假設角度以弧度的定點格式表示。要找到一個角度$\beta$的正弦函數和餘弦函數，可以在單位圓上找到對應角度的y座標和座標。利用CORDIC，會從以下的向量$v_0$開始：

$$v_0 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}.$$

在第一次的迭代時，向量逆時針轉45°，得到向量$v_1$。接著繼續的迭代，每一次的角度漸漸變小，旋轉方向可能順時針，也可能逆時針，直到得到想要的角度為止。每一次的角度為$\gamma_i = \arctan{(2^{-i})}$，其中$i = 0, 1, 2, \dots$。

若以更正式的方式表示，每一次的迭代就是一次旋轉，也就是將向量$v_i$乘以旋轉矩陣$R_{i}$：

$$v_{i+1} = R_i v_i.$$

旋轉矩陣為

$$R_i = \begin{bmatrix} <pre> \cos(\gamma_i) & -\sin(\gamma_i) \\ \sin(\gamma_i) & \cos(\gamma_i) </pre> \end{bmatrix}.$$

利用以下兩個三角恆等式：

$$\begin{align} <pre> \cos(\gamma_i) &= \frac{1}{\sqrt{1 + \tan^2(\gamma_i)}}, \\ \sin(\gamma_i) &= \frac{\tan(\gamma_i)}{\sqrt{1 + \tan^2(\gamma_i)}}, </pre> \end{align}$$

旋轉矩陣變成

$$R_i = \frac{1}{\sqrt{1 + \tan^2(\gamma_i)}} \begin{bmatrix} <pre> 1 & -\tan(\gamma_i) \\ \tan(\gamma_i) & 1 </pre> \end{bmatrix}.$$

旋轉向量$v_{i+1} = R_i v_i$就會變成下式

$$\begin{bmatrix} <pre> x_{i+1} \\ y_{i+1} </pre> \end{bmatrix} = \frac{1}{\sqrt{1 + \tan^2(\gamma_i)}} \begin{bmatrix} <pre> 1 & -\tan(\gamma_i) \\ \tan(\gamma_i) & 1 </pre> \end{bmatrix} \begin{bmatrix} <pre> x_i \\ y_i </pre> \end{bmatrix},$$

其中$x_i$和$y_i$是$v_i$的分量，若將角度$\gamma_i$限制在使$\tan(\gamma_i) = \pm 2^{-i}$的值，和tangent的乘法就以變成乘（或除）2的冪次，在數位電腦硬體中可以快速的用位元右移或左移來計算。因此上法會變成

$$\begin{bmatrix} <pre> x_{i+1} \\ y_{i+1} </pre> \end{bmatrix} = K_i \begin{bmatrix} <pre> 1 & -\sigma_i 2^{-i} \\ \sigma_i 2^{-i} & 1 </pre> \end{bmatrix} \begin{bmatrix} <pre> x_i \\ y_i </pre> \end{bmatrix},$$

其中

$$K_i = \frac{1}{\sqrt{1 + 2^{-2i}}},$$

且$\sigma_i$是用來判斷旋轉方向的。若角度$\gamma_i$為正，則$\sigma_i$為+1，否則則為−1。

所有的$K_i$因子可以在迭代過程中忽略，最後再一次乘以$K(n)$因子即可：

$$K(n) = \prod_{i=0}^{n-1} K_i = \prod_{i=0}^{n-1} \frac{1}{\sqrt{1 + 2^{-2i}}},$$

此數可以事先計算好存在表格中，若迭代次數是固定的，只需計算一個常數且儲存即可，甚至此修正也可以事先進行，將常數先乘以$v_0$，可以節省一次乘法。另外，可以注意到

$$K = \lim_{n \to \infty} K(n) \approx 0.6072529350088812561694$$

因此可以簡化演演算法的複雜度。有些應用會避免修正$K$，因此此演演算法本身會帶一個增益$A$：

$$A = \frac{1}{K} = \lim_{n \to \infty} \prod_{i=0}^{n-1} \sqrt{1 + 2^{-2i}} \approx 1.64676025812107.$$

在夠多次的迭代後，向量的角度會接近想要的角度$\beta$。以一般的應用來說，40次的迭代（n = 40）已可以有小數10位的精度。

唯一未解決的問題是判斷每一次迭代要順時針旋轉或逆時針旋轉（選擇$\sigma$的值）。這可以記錄每一次旋轉的角度，從還需要旋轉的角度中減去此一角度，會得到下一個還需要旋轉的角度$\beta$。若$\beta_{n+1}$為正，需要順時針旋轉，否則，就需要逆時針旋轉：

$$\beta_0 = \beta$$

$$\beta_{i+1} = \beta_i - \sigma_i \gamma_i, \quad \gamma_i = \arctan(2^{-i}).$$

$\gamma_n$的值需要事先計算且記錄下來。不過若是小角度，根據小角度近似，在定點數下可得$\arctan(\gamma_n) = \gamma_n$，因此可以節省儲存用的空間。

如以上所述，角度$\beta$的正弦函數為其y坐標，餘弦函數為其x坐標。

### 向量模式

上述旋轉模式的演演算法，是旋轉原來位在x軸上的單位向量。但此演演算法可以用來旋轉角度在−$\pi / 2$及$\pi / 2$之間的任意向量，旋轉的方向則依$\beta_i$的正負號決定。

向量模式下的演演算法和旋轉模式略有不同。其啟始向量的x坐標要為正值，y坐標則為任意值。持續轉動的目的是將向量旋轉到x軸（因此y座標為0）。每一步裡的旋轉方向會由y的值決定。$\beta_i$的最終值包括了總旋轉角度。x的最終值是原向量的大小，再乘以K。因此，可以看出向量模式可以進行直角坐標到極坐標的轉換。

## 實現

Java的Math類別中有`scalb(double x,int scale)`方法可以實現二進位的移位，C語言有ldexp函數，x86處理器有`fscale`浮點運算。

### 軟體範例（Python）

```python3
from math import atan2, sqrt, sin, cos, radians

ITERS = 16
theta_table = [atan2(1, 2**i) for i in range(ITERS)]

def compute_K(n):
    """
    Compute K(n) for n = ITERS. This could also be
    stored as an explicit constant if ITERS above is fixed.
    """
    k = 1.0
    for i in range(n):
        k *= 1 / sqrt(1 + 2 ** (-2 * i))
    return k

def CORDIC(alpha, n):
    K_n = compute_K(n)
    theta = 0.0
    x = 1.0
    y = 0.0
    P2i = 1  # This will be 2**(-i) in the loop below
    for arc_tangent in theta_table:
        sigma = +1 if theta < alpha else -1
        theta += sigma * arc_tangent
        x, y = x - sigma * y * P2i, sigma * P2i * x + y
        P2i /= 2
    return x * K_n, y * K_n

if __name__ == "__main__":
    # Print a table of computed sines and cosines, from -90° to +90°, in steps of 15°,
    # comparing against the available math routines.
    print("  x       sin(x)     diff. sine     cos(x)    diff. cosine ")
    for x in range(-90, 91, 15):
        cos_x, sin_x = CORDIC(radians(x), ITERS)
        print(
            f"{x:+05.1f}°  {sin_x:+.8f} ({sin_x-sin(radians(x)):+.8f}) {cos_x:+.8f} ({cos_x-cos(radians(x)):+.8f})"
        )
```

#### 輸出

```console
$ python cordic.py
  x       sin(x)     diff. sine     cos(x)    diff. cosine
-90.0°  -1.00000000 (+0.00000000) -0.00001759 (-0.00001759)
-75.0°  -0.96592181 (+0.00000402) +0.25883404 (+0.00001499)
-60.0°  -0.86601812 (+0.00000729) +0.50001262 (+0.00001262)
-45.0°  -0.70711776 (-0.00001098) +0.70709580 (-0.00001098)
-30.0°  -0.50001262 (-0.00001262) +0.86601812 (-0.00000729)
-15.0°  -0.25883404 (-0.00001499) +0.96592181 (-0.00000402)
+00.0°  +0.00001759 (+0.00001759) +1.00000000 (-0.00000000)
+15.0°  +0.25883404 (+0.00001499) +0.96592181 (-0.00000402)
+30.0°  +0.50001262 (+0.00001262) +0.86601812 (-0.00000729)
+45.0°  +0.70709580 (-0.00001098) +0.70711776 (+0.00001098)
+60.0°  +0.86601812 (-0.00000729) +0.50001262 (+0.00001262)
+75.0°  +0.96592181 (-0.00000402) +0.25883404 (+0.00001499)
+90.0°  +1.00000000 (-0.00000000) -0.00001759 (-0.00001759)
```

### 硬體範例

實現CORDIC需要的邏輯閘大約和乘法器相當，兩者都是用位元移位和加法所組合的。要選擇乘法器或是CORDIC會隨應用而定。若複數以其實部和虛部表示（直角座標），複數乘法會需要進行四次的乘法。但若複數以極座標表示，只要一個CORDIC即可處理，這更適合用在其乘積的量值不重要的應用（例如將向量和單位圓上的向量相乘的情形）。在數位下轉換器之類的通訊相關電路中，常會用到CORDIC。

## 雙重迭代CORDIC

在Vladimir Baykov的二篇文獻中，曾經提到用「雙重迭代」（double iterations）來實現反正弦函數、反餘弦函數、自然對數、指數函數以及雙曲函數。雙重迭代中的作法和傳統的CORDIC演演算法不同，傳統的CORDIC演演算法中，迭代變數會在每次迭代時加一。但在雙重迭代中，迭代變數會先重複一次，然後才加一。因此雙重迭代CORDIC演演算法的迭代變數會是$i = 0, 0, 1, 1, 2, 2\dots$，而傳統CORDIC演演算法的迭代變數是$i = 0, 1, 2\dots$。雙重迭代法保證方法的收斂，不過其引數的有效範圍會變化。

針對$R$進制數字系統中，通用的CORDIC收斂問題，可以證明針對正弦、餘弦及反正切函數，每一個i（i = 0或1到n，其中n是位數）進行$R - 1$次迭代就可以保證收斂。若是自然對數、指數、雙曲正弦、雙曲餘弦及雙曲反正切函數，每一個i需要進行$R$次迭代。若是反正弦函數及反餘弦函數，每一個i需要進行2 $R - 1$次的迭代。

若是雙曲反正弦及雙曲反餘弦函數，每一個i需要進行2R次的迭代。

## 相關演演算法

CORDIC屬於「移位及加法」（shift-and-add）的演演算法，就像Henry Briggs文獻提到的對數和指數演演算法一樣。BKM演演算法也是可以計算許多基本函數的移位及加法演演算法，是複數平面下對數和指數演演算法的推廣。BKM可以計算實數角度$x$（以弧度表示）的正弦和餘弦，其方式是計算$0+ix$的指數，也就是$\operatorname{cis}(x) = \cos(x) + i \sin(x)$。BKM演演算法比CORDIC要複雜一些，但其優點是不需考慮比例因子（K）。