mado: 三角函數

contributed by <jouae>

在 mado RP#27 中由ndsl7109256 貢獻改進三角函數查表的大小中，後續 ndsl7109256 採用五階固定點泰勒展開逼近正弦函數。

其中，遇到一個逼近問題「正切函數該怎麼計算？」

由於 mado 該專案希望最大限度地減少計算上使用到除法運算子 /，所以在計算正切函數時，避免直接以

$$\tan x={\displaystyle \frac{\sin x}{\cos x}}$$
計算。mado 中是希望減少除法運算子 / 的使用。

另一個計算三角函數可行的方案是採用座標旋轉數位計算法(Coordinate Rotation Digital Computer, CORDIC)。

該算法的優點在於僅需要加法和乘法就可以計算一三角函數的逼近值。

CORDIC 綜覽可以看 50 Years of CORDIC: Algorithms, Architectures, and Applications 一文

Coordinate Rotation Digital Computer

以浮點數實作部分

CORDIC 是 Jack E. Volder 在 1956 年提出的一份技術報告中的計算方法，並於 1959 年正式公開在 Proceedings of the Western Joint Computer Conference (WJCC) 上。

以下內容摘自原始論文 The CORDIC Computing Technique, Jack Volder, 1959

其精神在於對一向量做旋轉迭代。考慮一二維座標中的旋轉矩陣跟一二維向量

$[x_i,y_i]$，

$$\begin{array}{rl}\left[\begin{array}{c}{x}_{i+1}\\ y_{i+1}\end{array}\right]& =\left[\begin{array}{cc}\cos {\theta }_i& -\sin {\theta }_i\\ \sin {\theta }_i& \cos {\theta }_i\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]\\ & =\cos {\theta }_i\left[\begin{array}{cc}1& -\tan {\theta }_i\\ \tan {\theta }_i& 1\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]\end{array}$$

對於正整數

$i$ ，固定旋轉角度

${\theta }_i$ 使得

$\tan {\theta }_i=2^{-i}$ ，其中該角度

${\theta }_i$ 可以藉由反正切函數表示為

$${\theta }_i=\arctan (2^{-i})$$

則上述線性系統可以改寫成

$$\begin{array}{rl}\left[\begin{array}{c}{x}_{i+1}\\ y_{i+1}\end{array}\right]& =\cos (\arctan (2^{-i}))\left[\begin{array}{cc}1& -d_i{2}^{-i}\\ d_i{2}^{-i}& 1\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]\\ & =K_i\left[\begin{array}{cc}1& -d_i{2}^{-i}\\ d_i{2}^{-i}& 1\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]\end{array}$$
其中，

$K_i=\cos (\arctan (2^{-i}))$ 、

$d_i=1$ 當旋轉為順時針旋轉，

$d_i=-1$ 當旋轉為逆時針旋轉。藉由

$K_i$ 可以計算一縮放因子(scaling factor)

$K(n)$：

$$\begin{array}{rl}K(n)& =\prod _{i=0}^{n-1}{K}_i\\ & =\prod _{i=0}^{n-1}\cos (\arctan (2^{-i}))\\ & =\prod _{i=0}^{n-1}{\displaystyle \frac{1}{\sqrt{1+2^{-2i}}}}\end{array}$$

最後一個等式成立是基於

$\cos x=\frac{1}{\sqrt{1+{\tan }^{2}x}}$，此關係式從

$1+{\tan }^{2}x={\sec }^{2}x$ 推導。

$n$ 為一自設定的迭代次數。

具體來說把這些旋轉角度

${\theta }_i$ 加總起來則得到

$$\theta =\sum _{i=0}^{n-1}{d}_i{\theta }_i=\sum _{i=0}^{n-1}{d}_i\arctan (2^{-i})$$

其中

$\theta$ 為目標旋轉角度。

其中

$x_i$ 代表

$\cos {\theta }_i$ 的值，所求

$\cos \theta$ 為

$$x_{n-1}=\cos \theta =\cos \left(\sum _{i=0}^{n-1}{d}_i{\theta }_i\right)$$

而

$y_i$ 代表

$\sin {\theta }_i$ 的值，所求

$\sin \theta$ 為

$$y_{n-1}=\sin \theta =\sin \left(\sum _{i=0}^{n-1}{d}_i{\theta }_i\right)$$

為了確定

$d_i$ 為

$1$ 或

$-1$，令變數

$z_i$ 紀錄前

$i-1$ 次微旋轉角度總和，即

$$z_i=\sum _{k=0}^{i-1}{d}_k{\theta }_k,\text{for }0<i<n-1$$

則關係式可以表達為

$$z_{i+1}=z_i+d_i{\theta }_i$$

在令

$z=\theta$ 為目標旋轉角度，確認目標旋轉角度減去前

$i$ 次微旋轉角度總和的正負號即可得知

$d_i$

$$d_i=\{\begin{array}{ll}1,& \text{ if }z-z_i>0\\ -1,& \text{ if }z-z_i<0\end{array}$$

變數

$z_i$ 的作用很簡單，其實就是確定該向量在

$i$ 次旋轉後，下次該以逆時針(

$d_i=-1$)旋轉或是順時針旋轉(

$d_i=1$)朝目標角度旋轉該微旋轉角度

${\theta }_i$。

Max-Gulda 編寫的 Cordic-Math 以 C 語言實踐該算法。

收斂性分析

Expanding the range of convergence of the CORDIC algorithm, 1991

The Scaling Free Method

原始演算法的缺點之一在於需要計算縮放因子

$K_i$。

An array architecture for fast computation of discrete hartley transform, A. S. Dhar and S. Banerjee, 1991 一文中提出小角度近似三角函數計算旋轉矩陣的方法，並附上一誤差估計。

考慮正弦與餘弦的在點

$x_0+ϵ$ 泰勒展開

$$\begin{array}{rl}\sin (x_0+ϵ)& =\sin (x_0)+\frac{ϵ\cos (x_0)}{1!}+\frac{{ϵ}^2\sin (x_0)}{2!}-\frac{{ϵ}^3\cos (x_0)}{3!}+\dots \\ \cos (x_0+ϵ)& =\cos (x_0)+\frac{ϵ\sin (x_0)}{1!}-\frac{{ϵ}^2\cos (x_0)}{2!}+\frac{{ϵ}^3\sin (x_0)}{1!}+\dots \end{array}$$

當

$x_0=0$ 時，上述三角函數的低階逼近為

$$\begin{array}{rl}\sin (ϵ)& =ϵ-\frac{{ϵ}^2}{2!}\sin ({\xi }_s)\\ \cos (ϵ)& =1-\frac{{ϵ}^2}{2!}+\frac{{ϵ}^3}{3!}\sin ({\xi }_c)\end{array}$$

其中

$\frac{{ϵ}^2}{2!}\sin ({\xi }_s)$ 與

$\frac{{ϵ}^3}{3!}\sin ({\xi }_c)$ 為泰勒展開剩餘項，

${\xi }_s$ 和

${\xi }_c$ 皆介於

$0$ 至

$ϵ$ 之間。

簡易的說，用小角度

$ϵ$ 近似三角函數，如

$$\sin (ϵ)\approx ϵ,\cos (ϵ)\approx 1-\frac{{ϵ}^2}{2}$$

其誤差為

$$\begin{array}{rl}& |\sin (ϵ)-ϵ|=\frac{{ϵ}^2}{2}\sin ({\xi }_s)\le \frac{{ϵ}^2}{2}\\ & |\cos (ϵ)-(1-\frac{{ϵ}^2}{2})|=\frac{{ϵ}^3}{6}\sin ({\xi }_c)\le \frac{{ϵ}^3}{6}\end{array}$$

在此將小角度

$ϵ$ 皆以

$2^{-i}$ 取代，分析

$\frac{{ϵ}^3}{6}$ 則

$$\frac{2^{-3i}}{6}=2^{-3i}{2}^{-1}{3}^{-1}=2^{-(3i+1+{\log }_{2}3)}<2^{-(3i+2.585)}$$

最後不等式成立的原因為

$1+{\log }_{2}3\approx 2.58496250072$ 小於

$2.585$。表示餘弦函數小角度近似的誤差不超過

$2^{-(3i+2.585)}$：

$$|\cos (2^{-i})-(1-\frac{2^{-2i}}{2})|<2^{-(3i+2.585)}$$

又

$2^{-p}$ 表示右移

$p$ 位元運算，因為

$i$ 為一變動的數，想要讓精度保持在

$2^{-(3i+2.585)}$ 以內，令

$n$ 表最大迭代次數，即

$2^{-n}$ 為最小旋轉角度，則限制

$n$ 滿足：

$$2^{-n}\ge 2^{-(3i+2.585)}$$

上述不等式限制，是為了確保誤差項

$2^{-(3i+2.585)}$ 不超過最小旋轉角度

$2^{-n}$。

則

$$\lceil \frac{n-2.585}{3}\rceil \le i\le n$$

於是原旋轉矩陣可以近似表示為

$$\left[\begin{array}{cc}1-2^{-(2i+1)}& -2^{-i}\\ 2^{-i}& 1-2^{-(2i+1)}\end{array}\right]$$

所以最終旋轉後的座標表示為：

$$\begin{array}{rl}\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]& =\prod _{i=j}^{n-1}\left[\begin{array}{cc}1-2^{-(2i+1)}& -2^{-i}\\ 2^{-i}& 1-2^{-(2i+1)}\end{array}\right]\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]\end{array},\text{ for }j=\lceil \frac{p-2.585}{3}\rceil$$

The Look Up Table of the Angles

第

$i$ 微角度表示為

${\theta }_i=\arctan (2^{-i})$，則

$\tan {\theta }_i=2^{-i}$

$$\theta ={\theta }_L+{\theta }_S$$

https://www.researchgate.net/publication/382638647_Flexible_and_Cost-Effective_Spherical_to_Cartesian_Coordinate_Conversion_Using_3-D_CORDIC_Algorithm_on_FPGA

$i$	${\theta }_i$ (radix)	$\tan {\theta }_i$
$1$	$45$	$1$
$2$	$26.565$	$\frac{1}{2}$
$3$	$14.036$	$\frac{1}{4}$
$4$	$7.125$	$\frac{1}{8}$
$5$	$3.576$	$\frac{1}{16}$
$6$	$1.7876$	$\frac{1}{32}$
$7$	$0.8938$	$\frac{1}{64}$
$8$	$0.4469$	$\frac{1}{128}$
$9$	$0.2238$	$\frac{1}{256}$

以有限集合表示查找表

$S$ ，則，

$$S=\{\arctan 2^{-1},\arctan 2^{-2},\arctan 2^{-3},\dots ,\arctan 2^{-n}\}$$
其中，

$n$ 與經度

Hybird Method

1. 表格多大
2. 精度
3. 程式碼，相依性，一個 cycle 幾道指令

Reference

• https://github.com/Max-Gulda/Cordic-Math/blob/main/lib/cordicMath/src/cordic-math.c#L290-L332
• A Novel Scaling free Vectoring CORDIC and its FPGA Implementation
• CORDIC II: A New Improved CORDIC Algorithm
• http://zipcpu.com/dsp/2017/10/02/cordic-tb.html
• https://github.com/ZipCPU/cordic?tab=readme-ov-file
• Modified virtually scaling-free adaptive cordic rotator algorithm and architecture
• https://repositories.lib.utexas.edu/server/api/core/bitstreams/a2a8f7dc-4463-4f4b-a075-15a17b8e7efa/content
• An array architecture for fast computation of discrete hartley transform
• https://www.eng.auburn.edu/~agrawvd/COURSE/E7950_Fall09/TALK/cordic_presentation.pdf
• Expanding the range of convergence of the CORDIC algorithm