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 ¹ 算法之前，我一直以为计算器是利用泰勒公式去求解，最近才了解到原来在计算机中，CORDIC 算法远比泰勒公式高效。

下面我们就来了解下泰勒公式的不足之处和 CORDIC 算法是怎么做的。

泰勒公式的缺点

上一节我们使用泰勒公式去计算三角函数值时，需要做很多次乘法运算，而计算器中乘法运算是很昂贵的，其缺点如下所示：

展开过小则会导致展开精度不够，展开过大则会导致计算量过大；
幂运算需要使用乘法器，存在很多重复计算；
需要很多变量来存储中间值。

在之前的文章矩阵乘法的 Strassen 算法，就是通过减少乘法，多做加法，从而大大降低了运算量，那么我们可以用相同的思想来优化运算吗？

当然可以，让我们请出今天的主角：CORDIC 算法。

解析 CORDIC 算法

我们知道单位圆上一点 $P$ ，其坐标为： $(\cos \theta , \sin \theta )$ ，如下图所示：

Unit Cycle

如果我们接收一个旋转向量 $M_{in}$ 逆时针旋转 $\theta$ ，将点 $P(x_{in} , y_{in})$ 旋转到 $P'(x_{R} , y_{R})$ , 如下图所示：

CORDIC

很容易得到如下公式：

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} , i=0, 1,\dots , n

这样和 $\tan (\theta_{i})$ 乘法操作就变成了移位操作，我们知道计算机中移位操作是非常快的，就可以大大加快计算速度了。

但这里仍然有 $3$ 个问题需要解决：

对于任意角度，可以通过满足条件的角度累加来得到在数学上相同的结果吗？
每次旋转得到的结果仍然需要乘以 $\cos(\theta )$ ，这部分的计算成本如何？如何计算？
因为每次旋转角度 $\theta = \arctan(2^{-i})$ ，朝着目标角度进行旋转时，可能会出现没有超过目标角度的情况，也会存在超过目标角度的情况，这种情况如何解决呢？

CORDIC Expand

对于第一个问题，答案是否定的。可以从数学上证明只有 $\angle 45^{\circ}$ 的倍数角才可以得到完全一致的结果。但是在工程应用中，我们只需要满足一定精度即可，可以增加迭代次数无限逼近原始角度，如下所示提高 $n$ 值以无限逼近原始角度。

\theta_{d} = \sum_{i=0}^{n} \theta_{i} , \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}$