By Long Luo
任何一款科学计算器上都可以计算三角函数，三角函数应用于生活工作中的方方面面，但计算机是如何计算三角函数值的呢？
三角函数本质上是直角三角形的3条边的比例关系，计算并没有很直观，那么计算机是如何计算三角函数值的呢？
在微积分中我们学习过 泰勒公式 ，我们知道可以使用泰勒展开式来对某个值求得其近似值，例如：
\[ \sin \angle 18^{\circ} = \frac{\sqrt {5} - 1}{4} \approx 0.3090169943 \]
利用泰勒公式，取前 \(4\) 项：
\[ \sin x \approx x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} \]
代入可得：
\[ \sin \angle 18^{\circ} = \sin \frac{\pi}{10} = \frac{\pi}{10} - \frac{1}{6} (\frac{\pi}{10})^3 + \frac{1}{120} (\frac{\pi}{10})^5 - \frac{1}{5040} (\frac{\pi}{10})^7 \approx 0.30903399538 \]
可见取前 \(4\) 项时精度已经在 \(10^{-5}\) 之下，对于很多场合精度已经足够高了。
在没有了解 CORDIC(Coordinate Rotation Digital Computer) Algorithm 1 算法之前，我一直以为计算器是利用泰勒公式去求解，最近才了解到原来在计算机中，CORDIC 算法远比泰勒公式高效。
下面我们就来了解下泰勒公式的不足之处和 CORDIC 算法是怎么做的。
上一节我们使用泰勒公式去计算三角函数值时，需要做很多次乘法运算，而计算器中乘法运算是很昂贵的，其缺点如下所示：
在之前的文章 矩阵乘法的 Strassen 算法 ，就是通过减少乘法，多做加法，从而大大降低了运算量，那么我们可以用相同的思想来优化运算吗？
当然可以，让我们请出今天的主角：CORDIC 算法。
我们知道单位圆上一点 \(P\) ，其坐标为：\((\cos \theta , \sin \theta )\) ，如下图所示：
如果我们接收一个旋转向量 \(M_{in}\) 逆时针旋转 \(\theta\) ，将点 \(P(x_{in} , y_{in})\) 旋转到 \(P'(x_{R} , y_{R})\) , 如下图所示：
很容易得到如下公式：
\[ x_R = x_{in} \cos (\theta) - y_{in} \sin (\theta) \]
\[ y_R = x_{in} \sin (\theta) + y_{in} \cos (\theta) \]
实际上由 复数运算 ，我们知道复数乘法就是幅角相加，模长相乘。我们可以将上式写成下列矩阵运算形式：
\[ \begin{aligned} \begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix} = \begin{bmatrix} \cos (\theta) & -\sin (\theta) \\ \sin (\theta) & \cos (\theta ) \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix} \end{aligned} \]
但上式运算时，只是对向量 \(v_{in} = {\begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}}\) 进行了线性变换，乘以一个旋转向量 \(M_{in}\) ，得到了旋转后的结果：向量 \(v_{R} = {\begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix}}\) 。
但是上式仍然需要 \(4\) 次乘法和 \(2\) 次加减法操作， 复杂度没有任何降低，那怎么办呢？
当当…当！
通过上述分析，我们已经知道可以使用有限次旋转操作来避免复杂的乘法操作，我们修改矩阵运算公式，提取 \(\cos (\theta )\) ，则公式可以修改为：
\[ \begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix} = \cos (\theta) \begin{bmatrix} 1 & -\tan (\theta) \\ \tan (\theta) & 1 \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix} \]
如果我们选择合适的角度值 \(\theta_i\)，使得
\[ \tan (\theta_{i}) = 2^{-i} \quad i=0, 1,\dots , n \]
这样和 \(\tan (\theta_{i})\) 乘法操作就变成了移位操作，我们知道计算机中移位操作是非常快的，就可以大大加快计算速度了。
但这里仍然有 \(3\) 个问题需要解决：
对于第一个问题，答案是否定的。可以从数学上证明只有 \(\angle 45^{\circ}\) 的倍数角才可以得到完全一致的结果。但是在工程应用中，我们只需要满足一定精度即可，可以增加迭代次数无限逼近原始角度，如下所示提高 \(n\) 值以无限逼近原始角度。
\[ \theta_{d} = \sum_{i=0}^{n} \theta_{i} \quad \forall \tan (\theta_{i}) = 2^{-i} \]
对于第二个问题，我们先来个例子，以 \(57.535^{\circ}\) 为例来看看求解过程：
\[ 57.535^{\circ} = 45^{\circ} + 26.565^{\circ} - 14.03^{\circ} \]
那么第一次旋转：
\[ \begin{bmatrix} x_{0} \\ y_{0} \end{bmatrix} = \cos (45^{\circ}) \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix} \]
第二次旋转：
\[ \begin{bmatrix} x_{1} \\ y_{1} \end{bmatrix} = \cos (26.565^{\circ}) \begin{bmatrix} 1 & -2^{-1} \\ 2^{-1} & 1 \end{bmatrix} \begin{bmatrix} x_{0} \\ y_{0} \end{bmatrix} \]
第三次旋转：
\[ \begin{bmatrix} x_{2} \\ y_{2} \end{bmatrix} = \cos (-14.03^{\circ}) \begin{bmatrix} 1 & 2^{-2} \\ -2^{-2} & 1 \end{bmatrix} \begin{bmatrix} x_{1} \\ y_{1} \end{bmatrix} \]
综合可得：
\[ \begin{bmatrix} x_{2} \\ y_{2} \end{bmatrix} = \cos (45^{\circ}) \cos (26.565^{\circ}) \cos (-14.03^{\circ}) \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} 1 & -2^{-1} \\ 2^{-1} & 1 \end{bmatrix} \begin{bmatrix} 1 & 2^{-2} \\ -2^{-2} & 1 \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix} \]
因为 \(\tan (\theta_{i}) = 2^{-i}\) ，由三角公式可以计算出：
\[ \cos (\theta_{i}) = \frac {1}{\sqrt {1 + \tan ^{2}(\theta_{i})}} = \frac {1}{\sqrt {1 + 2^{-2i}}} \]
令 \(K_i = \cos (\theta_{i})\) ，则当进行 \(n\) 次迭代之后：
\[ K(n) = \prod _{i=0}^{n-1}K_{i} = \prod _{i=0}^{n-1}{\frac {1}{\sqrt {1 + 2^{-2i}}}} \]
当 \(\theta_{i}\) 越来越小时， \(\cos \theta\) 也越来越逼近 \(1\) ，当迭代次数 \(n \to \infty\) ， \(K(n)\) 极限存在，求解可得：
\[ K = \lim _{n \to \infty }K(n) \approx 0.6072529350088812561694 \]
由 \(K\) 我们实际上可以得到最终的向量 \(v_R\) 的模长极限为：
\[ A = \frac {1}{K} = \lim _{n \to \infty } \prod _{i=0}^{n - 1}{\sqrt {1 + 2^{-2i}}} \approx 1.64676025812107 \]
实际上当迭代次数为 \(6\) 时，可以计算出缩放比例 \(K\) ，就已经精确到 \(0.6072\) 了，如下所示：
\[ K \approx \cos (45^{\circ}) \cos (26.565^{\circ}) \times \dots \times \cos (0.895^{\circ}) = 0.6072 \]
实际上，任意角度只要迭代次数超过 \(6\) ，我们可以直接使用 \(K = 0.6072\) 这个值。
对于第三个问题，稍微有点复杂，我们在下一节继续讲解！
上一节遗留的问题是迭代旋转角度时，旋转角度不一定会落在目标角度内，我们需要引入一个角度误差，用来衡量旋转角度和目标角之间距离，如下所示：
\[ \theta_{error} = \theta_d - \sum_{i=0}^{n} \theta_{i} \]
当 \(\theta_{error} > 0\) 时，我们应该逆时针旋转，而 \(\theta_{error} < 0\) ，则顺时针旋转。根据精度需要，当 \(\left | \theta_{error} \right | \le \epsilon\) 即可退出迭代。
同时我们修改之前的公式，引入 \(\sigma_{i} \in \left \{ +1, -1 \right \}\) ，于是可以得到最终公式：
\[ x \left [ i+1 \right ] = x \left [ i \right ] - \sigma_{i} 2^{-i} y \left [ i \right ] \]
\[ y \left [ i+1 \right ] = y \left [ i \right ] + \sigma_{i} 2^{-i} x \left [ i \right ] \]
\[ z \left [ i+1 \right ] = z \left [ i \right ] - \sigma_{i} \tan^{-1} ( 2^{-i} ) \]
上面讲了这么多，来个实例吧，练习巩固下知识，看看自己是否真的懂了？
计算 \(\sin 70^{\circ}\) 和 \(\cos 70^{\circ}\) 的值。
从 \(x_{in}=1, y_{in}=0\) 开始，迭代 \(6\) 次结果如下：
| 第 \(i\) 次迭代 | \(\sigma_{i}\) | \(x \left[ i \right ]\) | \(y \left[ i \right ]\) | \(z \left[ i \right ]\) |
|---|---|---|---|---|
| - | - | \(1\) | \(0\) | \(70^{\circ}\) |
| \(0\) | \(1\) | \(1\) | \(1\) | \(25^{\circ}\) |
| \(1\) | \(1\) | \(0.5\) | \(1.5\) | \(-1.5651^{\circ}\) |
| \(2\) | \(-1\) | \(0.875\) | \(1.375\) | \(12.4711^{\circ}\) |
| \(3\) | \(1\) | \(0.7031\) | \(1.4844\) | \(5.3461^{\circ}\) |
| \(4\) | \(1\) | \(0.6103\) | \(1.5283\) | \(1.7698^{\circ}\) |
| \(5\) | \(1\) | \(0.5625\) | \(1.5474\) | \(-0.0201^{\circ}\) |
| \(6\) | \(-1\) | \(0.5867\) | \(1.5386\) | \(0.8751^{\circ}\) |
迭代到第 \(6\) 次时，角度误差已经小于 \(1^{\circ}\) 了， 通过表格可知：
\[ x_{R} = 0.6072 \times 0.5867 = 0.3562 \]
\[ y_{R} = 0.6072 \times 1.5386 = 0.9342 \]
通过计算器可知，\(\cos (70^{\circ}) = 0.34202\) ，\(\sin (70^{\circ}) = 0.93969\) ，误差已经在 \(\dfrac {1}{100}\) 之下了，实际应用中我们会迭代 \(16\) 次，误差会非常小。
通过上面分析，我们已经知道了 CORDIC 算法的原理，下面就开始编程吧！
用 JavaScript 写了一个在线互动版本，传送门 → ：
CORDIC 算法相比其他算法的优点，体现在以下几个方面：
几个月这篇文章发布之后，陆陆续续得到了不少读者的宝贵意见，有读者反馈说没有写 CORDIC 算法的缺点，也有读者反馈在他在嵌入式处理器STM32上还不如 泰勒展开式(Taylor Series) 快。
对于这 \(2\) 个问题，我又查阅了一些 CORDIC 算法2 的资料，可以回答下这 \(2\) 个问题。
CORDIC 算法的不足之处在于下面几点：
CORDIC 算法的收敛速度较慢，比如当输入值 \(\theta\) 很接近 \(0\) 或 \(\dfrac {\pi}{2}\) 时，这个时候就需要进行多轮迭代才能达到足够的精度。
当 \(\theta \approx 0\) 时，实际上只需要取前 \(2\) 项就有足够精度了，而且收敛速度快多了！
\[ \sin (x) = x - \frac {1}{3!}x^3 \]
\[ \cos (x) = 1 - \frac {1}{2!}x^2 \]
如上一小节，当输入值 \(\theta\) 是小角度或者大角度的情况时，CORDIC 算法的固定迭代次数可能也无法提供足够的精度。
CORDIC 算法在实际应用中需要考虑输入参数的限制 ，特别是幅值和角度范围，这是因为 CORDIC 算法是一种迭代算法，对于较大的角度，需要更多的迭代次数才能达到所需的精度，从而导致计算效率的下降。
CORDIC 算法适用于处理相对较小范围的角度，确切来说是在 \([−99.88^{\circ} , 99.88^{\circ}]\) 的范围内表现最佳，下面我们来证明下：
由 CORDIC 算法计算过程可知：
\[ \left | z[i] \right | \leq \sum_{j=i}^{\infty } \arctan 2^{-j} \]
那么可得：
\[ \theta_{max} = z[0]_{max} = \sum_{j=0}^{\infty } \arctan 2^{-j} \approx 1.7429(99.88^{\circ}) \quad \quad \sigma_j = 1, z[j] > 0 \]
假设 \(\theta < \theta_{max}\) ， 则有：
\[ \left | z[i] \right | \leq \arctan (2^{−(i−1)}) \]
因此
\[ \arctan (2^{−i}) \leq \sum_{j=i+1}^{\infty } \arctan (2^{−j}) \quad \forall i \]
所以 CORDIC算法的输入范围是：\([−99.88^{\circ} , 99.88^{\circ}]\) 。
对于第 \(2\) 个问题，我不太清楚那位读者是如何测试的！如之前所述，在某些情况下，泰勒展开式是优于 CORDIC 算法的！
CORDIC 算法是在硬件中实现的，是 FPGA 设计3 中的经典例题。
CORDIC 算法是一种高效计算三角函数值的方法。相比传统的泰勒展开式，它具有简单高效、低存储需求和迭代控制等优势。在需要计算三角函数值的应用中，CORDIC 算法更快速、更节省资源。
此内容由惯性聚合(RSS阅读器)自动聚合整理，仅供阅读参考。 原文来自 — 版权归原作者所有。