By Long Luo

任何一款科学核算器上都能够核算三角函数，三角函数运用于日子工作中的方方面面，但核算机是怎么核算三角函数值的呢？

三角函数本质上是直角三角形的3条边的份额关系，核算并没有很直观，那么核算机是怎么核算三角函数值的呢？

在微积分中咱们学习过泰勒公式，咱们知道能够运用泰勒展开式来对某个值求得其近似值，例如：

sin⁡∠18∘=5−14≈0.3090169943\sin \angle 18^{\circ} = \frac{\sqrt {5} – 1}{4} \approx 0.3090169943

运用泰勒公式，取前 $4$ 项：

sin⁡x≈x−x33!+x55!−x77!\sin x \approx x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!}

代入可得：

sin⁡∠18∘=sin⁡10=10−16(10)3+1120(10)5−15040(10)7≈0.30903399538\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⁡,sin⁡)(\cos \theta , \sin \theta )$ ，如下图所示：

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

很容易得到如下公式：

xR=xincos()−yinsin()x_R = x_{in} cos(\theta) – y_{in} sin(\theta)

yR=xinsin()+yincos()y_R = x_{in} sin(\theta) + y_{in} cos(\theta)

实践上由复数运算，咱们知道复数乘法便是幅角相加，模长相乘。咱们能够将上式写成下列矩阵运算形式：

[xRyR]=[cos⁡()−sin⁡()sin⁡()cos⁡()][xinyin]\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}

但上式运算时，仅仅对向量 $vin=[xinyin]v_{in} = {\begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}}$ 进行了线性变换，乘以一个旋转向量 $M_{in}$ ，得到了旋转后的成果：向量 $vR=[xRyR]v_{R} = {\begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix}}$ 。

但是上式依然需求 $4$ 次乘法和 $2$ 次加减法操作，复杂度没有任何下降，那怎样办呢？

当当…当！

经过上述剖析，咱们现已知道能够运用有限次旋转操作来避免复杂的乘法操作，咱们修改矩阵运算公式，提取 $cos⁡()\cos (\theta )$ ，则公式能够修改为：

[xRyR]=cos()[1−tan()tan()1][xinyin]\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}

如果咱们挑选合适的视点值 $i\theta_i$ ，使得

tan⁡(i)=2−i,i=0,1,…,n\tan (\theta_{i}) = 2^{-i} , i=0, 1,\dots , n

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

但这里依然有 $3$ 个问题需求处理：

关于恣意视点，能够经过满足条件的视点累加来得到在数学上相同的成果吗？
每次旋转得到的成果依然需求乘以 $cos⁡()\cos(\theta )$ ，这部分的核算成本怎么？怎么核算？
由于每次旋转视点 $=arctan⁡(2−i)\theta = \arctan(2^{-i})$ ，朝着方针视点进行旋转时，可能会出现没有超越方针视点的情况，也会存在超越方针视点的情况，这种情况怎么处理呢？

关于第一个问题，答案是否定的。能够从数学上证明只要 $∠45∘\angle 45^{\circ}$ 的倍数角才能够得到完全一致的成果。但是在工程运用中，咱们只需求满足一定精度即可，能够增加迭代次数无限迫临原始视点，如下所示进步 $n$ 值以无限迫临原始视点。

d=∑i=0ni,∀tan⁡(i)=2−i\theta_{d} = \sum_{i=0}^{n} \theta_{i} , \forall \tan(\theta_{i}) = 2^{-i}

关于第二个问题，咱们先来个例子，以 $57.535∘57.535^{\circ}$ 为例来看看求解进程：

57.535∘=45∘+26.565∘−14.03∘57.535^{\circ} = 45^{\circ}+26.565^{\circ}-14.03^{\circ}

那么第一次旋转：

[x0y0]=cos(45∘)[1−111][xinyin]\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}

第二次旋转：

[x1y1]=cos(26.565∘)[1−2−12−11][x0y0]\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}

第三次旋转：

[x2y2]=cos(−14.03∘)[12−2−2−21][x1y1]\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}

归纳可得：

[x2y2]=cos(45∘)cos(26.565∘)cos(−14.03∘)[1−111][1−2−12−11][12−2−2−21][xinyin]\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⁡(i)=2−i\tan (\theta_{i}) = 2^{-i}$ ，由三角公式能够核算出：

cos⁡(i)=11+tan⁡2(i)=11+2−2i\cos(\theta_{i}) = \frac {1}{\sqrt {1 + \tan ^{2}(\theta_{i})}} = \frac {1}{\sqrt {1 + 2^{-2i}}}

令 $Ki=cos⁡(i)K_i = \cos(\theta_{i})$ ，则当进行 $n$ 次迭代之后：

\prod _{i=0}^{n-1}K_{i} = \prod _{i=0}^{n-1}{\frac {1}{\sqrt {1 + 2^{-2i}}}}

当 $i\theta_{i}$ 越来越小时， $cos⁡\cos \theta$ 也越来越迫临 $1$ ，当迭代次数 $\to \infty$ ， $K (n)$ 极限存在，求解可得：

\lim _{n \to \infty }K(n) \approx 0.6072529350088812561694

由 $K$ 咱们实践上能够得到最终的向量 $v_R$ 的模长极限为：

{\frac {1}{K}} = \lim _{n \to \infty } \prod _{i=0}^{n – 1}{\sqrt {1 + 2^{-2i}}} \approx 1.64676025812107

实践受骗迭代次数为 $6$ 时，能够核算出缩放份额 $K$ ，就现已准确到 $0.6072$ 了，如下所示：

\approx cos(45^{\circ}) cos(26.565^{\circ}) \times \dots \times cos(0.895^{\circ}) = 0.6072

实践上，恣意视点只要迭代次数超越 $6$ ，咱们能够直接运用 $K = 0.6072$ 这个值。

关于第三个问题，略微有点复杂，咱们在下一节持续解说！

视点累加

上一节留传的问题是迭代旋转视点时，旋转视点纷歧定会落在方针视点内，咱们需求引进一个视点差错，用来衡量旋转视点和方针角之间间隔，如下所示：

error=d−∑i=0ni\theta_{error} = \theta_d – \sum_{i=0}^{n} \theta_{i}

当 $error>0\theta_{error} > 0$ 时，咱们应该逆时针旋转，而 $error<0\theta_{error} < 0$ ，则顺时针旋转。根据精度需求，当 $∣error∣≤\left | \theta_{error} \right | \le \epsilon$ 即可退出迭代。

同时咱们修改之前的公式，引进 $i∈{+1,−1}\sigma_{i} \in \left \{ +1, -1 \right \}$ ，于是能够得到最终公式：

\left [ i+1 \right ] = x \left [ i \right ] – \sigma_{i} 2^{-i} y \left [ i \right ]

\left [ i+1 \right ] = y \left [ i \right ] + \sigma_{i} 2^{-i} x \left [ i \right ]

\left [ i+1 \right ] = z \left [ i \right ] – \sigma_{i} tan^{-1} ( 2^{-i} )

举个例子

上面讲了这么多，来个实例吧，练习巩固下知识，看看自己是否真的懂了？

核算 $sin⁡70∘\sin 70^{\circ}$ 和 $cos⁡70∘\cos 70^{\circ}$ 的值。

从 $x_{in}=1, y_{in} = 0$ 开端，迭代 $6$ 次成果如下：

第 $i$ 次迭代	$i\sigma_{i}$	$\left[ i \right ]$	$\left[ i \right ]$	$\left[ i \right ]$
–	–	$1$	$0$	$70∘70^{\circ}$
$0$	$1$	$1$	$1$	$25∘25^{\circ}$
$1$	$1$	$0.5$	$1.5$	$−1.5651∘-1.5651^{\circ}$
$2$	$- 1$	$0.875$	$1.375$	$12.4711∘12.4711^{\circ}$
$3$	$1$	$0.7031$	$1.4844$	$5.3461∘5.3461^{\circ}$
$4$	$1$	$0.6103$	$1.5283$	$1.7698∘1.7698^{\circ}$
$5$	$1$	$0.5625$	$1.5474$	$−0.0201∘-0.0201^{\circ}$
$6$	$- 1$	$0.5867$	$1.5386$	$0.8751∘0.8751^{\circ}$

迭代到第 $6$ 次时，视点差错现已小于 $1∘1^{\circ}$ 了，经过表格可知：

xR=0.60720.5867=0.3562x_{R} = 0.6072 \times 0.5867 = 0.3562

yR=0.60721.5386=0.9342y_{R} = 0.6072 \times 1.5386 = 0.9342

经过核算器可知， $cos⁡(70∘)=0.34202\cos(70^{\circ}) = 0.34202$ ， $sin⁡(70∘)=0.93969\sin(70^{\circ}) = 0.93969$ ，差错现已在 $1100\frac {1}{100}$ 之下了，实践运用中咱们会迭代 $16$ 次，差错会非常小。

在线 CORDIC 算法Demo

经过上面剖析，咱们现已知道了 CORDIC 算法的原理，下面就开端编程吧！

用 JavaScript 写了一个在线互动版别，传送门 → ：

能够调整不同迭代次数，和系统库函数 $Math.cos\textit{Math.cos}$ ， $Math.sin\textit{Math.sin}$ 进行比较：

能够检查每次迭代的成果，把握 Cordic 算法迭代原理：

CORDIC 算法的长处

CORDIC 算法相比其他算法的长处，体现在以下几个方面：

简化运算：CORDIC 算法首要运用位移、加法和减法等简单的运算，避免了复杂的乘法操作，从而进步了核算速度；
并行核算：CORDIC 算法的迭代操作是相互独立的，能够进行并行核算，运用现代核算机的多核优势，进一步提高核算功率；
硬件优化：CORDIC 算法适用于硬件实现，能够经过专用硬件电路（如FPGA）进行加快，极大地进步核算速度；
低存储需求：CORDIC 算法只需存储一小组预先核算好的旋转视点，节省了存储空间；
迭代操控：经过操控迭代次数，能够平衡核算精度和核算速度，根据需求进行调整。

总结

CORDIC 算法是一种高效核算三角函数值的办法。相比传统的泰勒展开式，它具有简单高效、低存储需求和迭代操控等优势。在需求核算三角函数值的运用中，CORDIC 算法更快速、更节省资源。

博客网站原文链接：www.longluo.me/blog/2023/0…

参考文献

CORDIC ↩

CORDIC算法：一种高效计算三角函数值的方法

泰勒公式的缺点

解析 CORDIC 算法

视点累加

举个例子

在线 CORDIC 算法Demo

CORDIC 算法的长处

总结

参考文献

作者信息

CORDIC算法：一种高效计算三角函数值的方法

泰勒公式的缺点

解析 CORDIC 算法

视点累加

举个例子

在线 CORDIC 算法Demo

CORDIC 算法的长处

总结

参考文献

Footnotes

相关文章

Matlab常用图像处理命令108例（六）

从零开始搭建个人网站

美团热更新方案Robust介绍

分布式存储系统 Ceph 介绍与环境部署

作者信息