By Long Luo
任何一款科学核算器上都能够核算三角函数,三角函数运用于日子工作中的方方面面,但核算机是怎么核算三角函数值的呢?
三角函数本质上是直角三角形的3条边的份额关系,核算并没有很直观,那么核算机是怎么核算三角函数值的呢?
在微积分中咱们学习过 泰勒公式 ,咱们知道能够运用泰勒展开式来对某个值求得其近似值,例如:
sin∠18∘=5−14≈0.3090169943\sin \angle 18^{\circ} = \frac{\sqrt {5} – 1}{4} \approx 0.3090169943sin∠18∘=45−1≈0.3090169943
运用泰勒公式,取前 444 项:
sinx≈x−x33!+x55!−x77!\sin x \approx x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!}sinx≈x−3!x3+5!x5−7!x7
代入可得:
sin∠18∘=sin10=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.30903399538sin∠18∘=sin10=10−61(10)3+1201(10)5−50401(10)7≈0.30903399538
可见取前 444 项时精度现已在 10−510^{-5}10−5 之下,关于许多场合精度现已足够高了。
在没有了解 CORDIC(Coordinate Rotation Digital Computer) Algorithm 1 算法之前,我一向以为核算器是运用泰勒公式去求解,最近才了解到原来在核算机中,CORDIC 算法远比泰勒公式高效。
下面咱们就来了解下泰勒公式的不足之处和 CORDIC 算法是怎样做的。
泰勒公式的缺点
上一节咱们运用泰勒公式去核算三角函数值时,需求做许多次乘法运算,而核算器中乘法运算是很昂贵的,其缺点如下所示:
- 展开过小则会导致展开精度不够,展开过大则会导致核算量过大;
- 幂运算需求运用乘法器,存在许多重复核算;
- 需求许多变量来存储中心值。
在之前的文章 矩阵乘法的 Strassen 算法 ,便是经过削减乘法,多做加法,从而大大下降了运算量,那么咱们能够用相同的思想来优化运算吗?
当然能够,让咱们请出今日的主角:CORDIC 算法。
解析 CORDIC 算法
咱们知道单位圆上一点 PPP ,其坐标为:(cos,sin)(\cos \theta , \sin \theta )(cos,sin) ,如下图所示:
如果咱们接收一个旋转向量 MinM_{in}Min 逆时针旋转 \theta ,将点 P(xin,yin)P(x_{in} , y_{in})P(xin,yin) 旋转到 P′(xR,yR)P'(x_{R} , y_{R})P′(xR,yR) , 如下图所示:
很容易得到如下公式:
xR=xincos()−yinsin()x_R = x_{in} cos(\theta) – y_{in} sin(\theta)xR=xincos()−yinsin()
yR=xinsin()+yincos()y_R = x_{in} sin(\theta) + y_{in} cos(\theta)yR=xinsin()+yincos()
实践上由 复数运算 ,咱们知道复数乘法便是幅角相加,模长相乘。咱们能够将上式写成下列矩阵运算形式:
[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}[xRyR]=[cos()sin()−sin()cos()][xinyin]
但上式运算时,仅仅对向量 vin=[xinyin]v_{in} = {\begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}}vin=[xinyin] 进行了线性变换,乘以一个旋转向量 MinM_{in}Min ,得到了旋转后的成果:向量 vR=[xRyR]v_{R} = {\begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix}}vR=[xRyR] 。
但是上式依然需求 444 次乘法和 222 次加减法操作, 复杂度没有任何下降,那怎样办呢?
当当…当!
经过上述剖析,咱们现已知道能够运用有限次旋转操作来避免复杂的乘法操作,咱们修改矩阵运算公式,提取 cos()\cos (\theta )cos() ,则公式能够修改为:
[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}[xRyR]=cos()[1tan()−tan()1][xinyin]
如果咱们挑选合适的视点值 i\theta_ii,使得
tan(i)=2−i,i=0,1,…,n\tan (\theta_{i}) = 2^{-i} , i=0, 1,\dots , ntan(i)=2−i,i=0,1,…,n
这样和 tan(i)\tan (\theta_{i})tan(i) 乘法操作就变成了移位操作,咱们知道核算机中移位操作是非常快的,就能够大大加快核算速度了。
但这里依然有 333 个问题需求处理:
- 关于恣意视点 ,能够经过满足条件的视点累加来得到在数学上相同的成果吗?
- 每次旋转得到的成果依然需求乘以 cos()\cos(\theta )cos() ,这部分的核算成本怎么?怎么核算?
- 由于每次旋转视点 =arctan(2−i)\theta = \arctan(2^{-i})=arctan(2−i) ,朝着方针视点进行旋转时,可能会出现没有超越方针视点的情况,也会存在超越方针视点的情况,这种情况怎么处理呢?
关于第一个问题,答案是否定的。能够从数学上证明只要 ∠45∘\angle 45^{\circ}∠45∘ 的倍数角才能够得到完全一致的成果。但是在工程运用中,咱们只需求满足一定精度即可,能够增加迭代次数无限迫临原始视点,如下所示进步 nnn 值以无限迫临原始视点。
d=∑i=0ni,∀tan(i)=2−i\theta_{d} = \sum_{i=0}^{n} \theta_{i} , \forall \tan(\theta_{i}) = 2^{-i}d=i=0∑ni,∀tan(i)=2−i
关于第二个问题,咱们先来个例子,以 57.535∘57.535^{\circ}57.535∘ 为例来看看求解进程:
57.535∘=45∘+26.565∘−14.03∘57.535^{\circ} = 45^{\circ}+26.565^{\circ}-14.03^{\circ}57.535∘=45∘+26.565∘−14.03∘
那么第一次旋转:
[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}[x0y0]=cos(45∘)[11−11][xinyin]
第二次旋转:
[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}[x1y1]=cos(26.565∘)[12−1−2−11][x0y0]
第三次旋转:
[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(−14.03∘)[1−2−22−21][x1y1]
归纳可得:
[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}[x2y2]=cos(45∘)cos(26.565∘)cos(−14.03∘)[11−11][12−1−2−11][1−2−22−21][xinyin]
由于 tan(i)=2−i\tan (\theta_{i}) = 2^{-i}tan(i)=2−i ,由三角公式能够核算出:
cos(i)=11+tan2(i)=11+2−2i\cos(\theta_{i}) = \frac {1}{\sqrt {1 + \tan ^{2}(\theta_{i})}} = \frac {1}{\sqrt {1 + 2^{-2i}}}cos(i)=1+tan2(i)1=1+2−2i1
令 Ki=cos(i)K_i = \cos(\theta_{i})Ki=cos(i) ,则当进行 nnn 次迭代之后:
K(n)=∏i=0n−1Ki=∏i=0n−111+2−2i K(n) = \prod _{i=0}^{n-1}K_{i} = \prod _{i=0}^{n-1}{\frac {1}{\sqrt {1 + 2^{-2i}}}}K(n)=i=0∏n−1Ki=i=0∏n−11+2−2i1
当 i\theta_{i}i 越来越小时, cos\cos \thetacos 也越来越迫临 111 ,当迭代次数 n→∞n \to \inftyn→∞ , K(n)K(n)K(n) 极限存在,求解可得:
K=limn→∞K(n)≈0.6072529350088812561694K = \lim _{n \to \infty }K(n) \approx 0.6072529350088812561694K=n→∞limK(n)≈0.6072529350088812561694
由 KKK 咱们实践上能够得到最终的向量 vRv_RvR 的模长极限为:
A=1K=limn→∞∏i=0n−11+2−2i≈1.64676025812107A = {\frac {1}{K}} = \lim _{n \to \infty } \prod _{i=0}^{n – 1}{\sqrt {1 + 2^{-2i}}} \approx 1.64676025812107A=K1=n→∞limi=0∏n−11+2−2i≈1.64676025812107
实践受骗迭代次数为 666 时,能够核算出缩放份额 KKK ,就现已准确到 0.60720.60720.6072 了,如下所示:
K≈cos(45∘)cos(26.565∘)⋯cos(0.895∘)=0.6072K \approx cos(45^{\circ}) cos(26.565^{\circ}) \times \dots \times cos(0.895^{\circ}) = 0.6072K≈cos(45∘)cos(26.565∘)⋯cos(0.895∘)=0.6072
实践上,恣意视点只要迭代次数超越 666 ,咱们能够直接运用 K=0.6072K = 0.6072K=0.6072 这个值。
关于第三个问题,略微有点复杂,咱们在下一节持续解说!
视点累加
上一节留传的问题是迭代旋转视点时,旋转视点纷歧定会落在方针视点内,咱们需求引进一个视点差错,用来衡量旋转视点和方针角之间间隔,如下所示:
error=d−∑i=0ni\theta_{error} = \theta_d – \sum_{i=0}^{n} \theta_{i}error=d−i=0∑ni
当 error>0\theta_{error} > 0error>0 时,咱们应该逆时针旋转,而 error<0\theta_{error} < 0error<0 ,则顺时针旋转。根据精度需求,当 ∣error∣≤\left | \theta_{error} \right | \le \epsilon∣error∣≤ 即可退出迭代。
同时咱们修改之前的公式,引进 i∈{+1,−1}\sigma_{i} \in \left \{ +1, -1 \right \}i∈{+1,−1} ,于是能够得到最终公式:
x[i+1]=x[i]−i2−iy[i]x \left [ i+1 \right ] = x \left [ i \right ] – \sigma_{i} 2^{-i} y \left [ i \right ]x[i+1]=x[i]−i2−iy[i]
y[i+1]=y[i]+i2−ix[i]y \left [ i+1 \right ] = y \left [ i \right ] + \sigma_{i} 2^{-i} x \left [ i \right ]y[i+1]=y[i]+i2−ix[i]
z[i+1]=z[i]−itan−1(2−i)z \left [ i+1 \right ] = z \left [ i \right ] – \sigma_{i} tan^{-1} ( 2^{-i} )z[i+1]=z[i]−itan−1(2−i)
举个例子
上面讲了这么多,来个实例吧,练习巩固下知识,看看自己是否真的懂了?
核算 sin70∘\sin 70^{\circ}sin70∘ 和 cos70∘\cos 70^{\circ}cos70∘ 的值。
从 xin=1,yin=0x_{in}=1, y_{in} = 0xin=1,yin=0 开端,迭代 666 次成果如下:
第 iii 次迭代 |
i\sigma_{i}i |
x[i]x \left[ i \right ]x[i] |
y[i]y \left[ i \right ]y[i] |
z[i]z \left[ i \right ]z[i] |
– |
– |
111 |
000 |
70∘70^{\circ}70∘ |
000 |
111 |
111 |
111 |
25∘25^{\circ}25∘ |
111 |
111 |
0.50.50.5 |
1.51.51.5 |
−1.5651∘-1.5651^{\circ}−1.5651∘ |
222 |
−1-1−1 |
0.8750.8750.875 |
1.3751.3751.375 |
12.4711∘12.4711^{\circ}12.4711∘ |
333 |
111 |
0.70310.70310.7031 |
1.48441.48441.4844 |
5.3461∘5.3461^{\circ}5.3461∘ |
444 |
111 |
0.61030.61030.6103 |
1.52831.52831.5283 |
1.7698∘1.7698^{\circ}1.7698∘ |
555 |
111 |
0.56250.56250.5625 |
1.54741.54741.5474 |
−0.0201∘-0.0201^{\circ}−0.0201∘ |
666 |
−1-1−1 |
0.58670.58670.5867 |
1.53861.53861.5386 |
0.8751∘0.8751^{\circ}0.8751∘ |
迭代到第 666 次时,视点差错现已小于 1∘1^{\circ}1∘ 了, 经过表格可知:
xR=0.60720.5867=0.3562x_{R} = 0.6072 \times 0.5867 = 0.3562xR=0.60720.5867=0.3562
yR=0.60721.5386=0.9342y_{R} = 0.6072 \times 1.5386 = 0.9342yR=0.60721.5386=0.9342
经过核算器可知,cos(70∘)=0.34202\cos(70^{\circ}) = 0.34202cos(70∘)=0.34202 ,sin(70∘)=0.93969\sin(70^{\circ}) = 0.93969sin(70∘)=0.93969 ,差错现已在 1100\frac {1}{100}1001 之下了,实践运用中咱们会迭代 161616 次,差错会非常小。
在线 CORDIC 算法Demo
经过上面剖析,咱们现已知道了 CORDIC 算法的原理,下面就开端编程吧!
用 JavaScript 写了一个在线互动版别,传送门 → :
- 能够调整不同迭代次数,和系统库函数 Math.cos\textit{Math.cos}Math.cos ,Math.sin\textit{Math.sin}Math.sin 进行比较:
- 能够检查每次迭代的成果,把握 Cordic 算法迭代原理:
CORDIC 算法的长处
CORDIC 算法相比其他算法的长处,体现在以下几个方面:
- 简化运算:CORDIC 算法首要运用位移、加法和减法等简单的运算,避免了复杂的乘法操作,从而进步了核算速度;
- 并行核算:CORDIC 算法的迭代操作是相互独立的,能够进行并行核算,运用现代核算机的多核优势,进一步提高核算功率;
- 硬件优化:CORDIC 算法适用于硬件实现,能够经过专用硬件电路(如FPGA)进行加快,极大地进步核算速度;
- 低存储需求:CORDIC 算法只需存储一小组预先核算好的旋转视点,节省了存储空间;
- 迭代操控:经过操控迭代次数,能够平衡核算精度和核算速度,根据需求进行调整。
总结
CORDIC 算法是一种高效核算三角函数值的办法。相比传统的泰勒展开式,它具有简单高效、低存储需求和迭代操控等优势。在需求核算三角函数值的运用中,CORDIC 算法更快速、更节省资源。
博客网站原文链接:www.longluo.me/blog/2023/0…
参考文献
-
CORDIC ↩