前言
CORDIC 算法比我想的还要难一点,也分多篇来写。
CORDIC 算法
基本思想
CORDIC 算法的基本原理是旋转向量逼近待求角 $\theta$ ,再根据直角三角形上三角函数的定义,将向量分量作为 sin/cos 值返回。
向量旋转
将平面上的向量 $\vec{V}=(x,y)$ 旋转 $\theta$ 角,公式为
$$ \vec{V_r} = \begin{cases} x_r = x \cdot \cos\theta - y \cdot \sin\theta \\ y_r = y \cdot \cos\theta + x \cdot \sin\theta \end{cases} $$
将 $\cos\theta$ 提取公因式,由于 $\tan\theta = \frac{\sin\theta}{\cos\theta}$,得到
$$ \vec{V_r} = \begin{cases} x_r = \cos\theta (x - y \cdot \tan\theta) \\ y_r = \cos\theta (y + x \cdot \tan\theta) \end{cases} $$
迭代逼近
CORDIC 是一种迭代逼近算法,它将原始的 $\theta$ 角拆分成 $N$ 次旋转操作,每次旋转一个细微角度来逼近:$\theta = \sum\limits_{i=0}^{N-1}{d_i \cdot a_i}$
$d_i \cdot a_i$ 表示的是第 $i$ 次旋转的弧度值,$d_i$ 是旋转方向,正为逆时针,负为顺时针。
我们就可以写出迭代公式:
$$ \vec{V_{i+1}} = \begin{cases} x_{i+1} = \cos(d_i \cdot a_i) (x_i - y_i \cdot \tan(d_i \cdot a_i)) \\ y_{i+1} = \cos(d_i \cdot a_i) (y_i + x_i \cdot \tan(d_i \cdot a_i)) \end{cases} $$
又因为 $d_i$ 实际上取值只有正负号,$\tan$ 是奇函数,$\tan(-x)=-\tan(x)$,化简得到
$$ \vec{V_{i+1}} = \begin{cases} x_{i+1} = \cos(d_i \cdot a_i) (x_i - y_i \cdot d_i \cdot \tan(a_i)) \\ y_{i+1} = \cos(d_i \cdot a_i) (y_i + x_i \cdot d_i \cdot \tan(a_i)) \end{cases} $$
我们再考虑:$\cos(d_i \cdot a_i)$
因为 $\cos$ 是偶函数,关于 y 轴对称,正负号不影响取值。因此,我们可以从 $\cos$ 里移除 $d_i$ 系数。
$$ \vec{V_{i+1}} = \begin{cases} x_{i+1} = \cos(a_i) (x_i - y_i \cdot d_i \cdot \tan(a_i)) \\ y_{i+1} = \cos(a_i) (y_i + x_i \cdot d_i \cdot \tan(a_i)) \end{cases} $$
选择特殊弧度值
在每一轮迭代中,选择特殊的旋转弧度,让 $\tan(a_i)$ 的取值为 $2^{-i}$ ,也就是设 $a_i = \arctan(2^{-i})$
从而让迭代公式变为
$$ \vec{V_{i+1}} = \begin{cases} x_{i+1} = \cos(a_i) \cdot (x_i - y_i \cdot d_i \cdot 2^{-i}) \\ y_{i+1} = \cos(a_i) \cdot (y_i + x_i \cdot d_i \cdot 2^{-i}) \end{cases} $$
提取 cos(ai) 系数
我们以计算 $x_{i+2}$ 为例,观察展开过程:
$$ \begin{aligned} x_{i+2} &= \cos(a_{i+1}) \cdot (x_{i+1} - y_{i+1} \cdot d_{i+1} \cdot 2^{-i+1}) \\ &= \cos(a_{i+1}) \cdot [ \cos(a_i) \cdot (x_i - y_i \cdot d_i \cdot 2^{-i}) - \cos(a_i) \cdot (y_i + x_i \cdot d_i \cdot 2^{-i}) \cdot d_{i+1} \cdot 2^{-i+1} ] \\ &= \cos(a_{i+1}) \cdot \cos(a_i) \cdot [ (x_i - y_i \cdot d_i \cdot 2^{-i}) - (y_i + x_i \cdot d_i \cdot 2^{-i}) \cdot d_{i+1} \cdot 2^{-i+1} ] \\ \end{aligned} $$
每次展开就会多一个 $\cos(a_i)$,右边括号里则多一个加减法和位移操作。
从 $x_i,y_i$ 的定义里不难发现发现,每个分量都会展开得到一个 $\cos(a)$ 的因子。
所以进一步展开,我们依然可以得到新的公因式。
此处可以用矩阵乘法的方式更形式化地证明可以提取出 $\prod^{N-1}_{i=0}\cos(a_i)$ 。
还是问 AI 要证明过程吧。
我们可以看到 $cos(a_i)$ 这个公因式是可以被提取出来,迭代 $N$ 轮就是叠乘 $N$ 次,记为:
$$ \begin{aligned} K &= \cos(a_0) \cdot \cos(a_1) \cdot \cos(a_2) \cdot ... \cdot \cos(a_i) \\ &= \prod_{i=0}^{N-1}\cos(a_i) \end{aligned} $$
所以 $\cos(a_i)$ 我们都可以提取到迭代最后进行计算,N轮迭代得到:
$$ \vec{V_{r}} = \begin{cases} x_{r} = K \cdot x_N \\ y_{r} = K \cdot y_N \end{cases} $$
其中 $x_N, y_N$ 都是没有 $\cos(a_i)$ 系数的迭代结果,提取公因式后的新的迭代公式为:
$$ \begin{cases} x_{i+1} = x_i - y_i \cdot d_i \cdot 2^{-i} \\ y_{i+1} = y_i + x_i \cdot d_i \cdot 2^{-i} \end{cases} $$
求系数 K
因为我们预设了 $2^{-i}=\tan(a_i)$ ,已知 $a_i = \arctan(2^{-i})$,代入可得:
$$ \begin{aligned} K &= \prod_{i=0}^{N-1} \cos(a_i) \\ &= \prod_{i=0}^{N-1} \cos(\arctan(2^{-i})) \end{aligned} $$
程序中我们一般迭代的是一个固定次数 $N$,可以预先求出 $K$ 的值。
实际上,当 $N=\infty$ 时,K 会收敛到一个极限值。具体的证明过程可以问 AI 。
所以我们其实也可以将 $K$ 直接取这个极限值,无论迭代多少次都乘这个固定系数。
我们使用 Python 计算 $K$ 的近似数值解:
import math
def calculate_K(N):
K = 1.0 # 初始化 K 为 1.0,因为是乘积
for i in range(N): # 循环从 0 到 N (不包括 N)
term = math.cos(math.atan(2**(-i)))
K *= term
return K
N_value = 5
result_K = calculate_K(N_value)
print(f"当 N = {N_value} 时,K 的值为: {result_K:.09f}")
N_value_large = 10
result_K_large = calculate_K(N_value_large)
print(f"当 N = {N_value_large} 时,K 的值为: {result_K_large:.09f}")
N_value_large = 30
result_K_large = calculate_K(N_value_large)
print(f"当 N = {N_value_large} 时,K 的值为: {result_K_large:.09f}")
当 N = 5 时,K 的值为: 0.607351770
当 N = 10 时,K 的值为: 0.607253032
当 N = 30 时,K 的值为: 0.607252935
我们取 $K \approx 0.607252935$ 来计算。
待续
下一部分会尝试把 CORDIC 算法推广到定点数,用 go 实现常用定点数三角函数 sin
/cos
/tan
。