基于CORDIC算法实现sin()函数（C语言）

CORDIC算法简介：

CORDIC（Coordinate Rotation Digital Computer）坐标旋转数学计算是数学与计算机技术交叉之后所产生的一中机器学习算法，其广泛应用于解决计算机的数学计算问题中。发展到现在，CORDIC算法及其扩展算法大致有三种计算模式：圆周旋转模式、线性旋转模式和双曲线旋转模式，分别用来实现不同的数学运算。

简单地说，CORDIC利用近似逼近的思想，将计算机中三角函数、开根号、求对数等复杂运算，转化为简单的加减和移位操作。

（1）算法原理

通过将原坐标轴转换为变换后的坐标轴后，计算转角

（2）算法理论分析

假设在原坐标系中，A的坐标为(X,Y),记在转换后的坐标中A的坐标为(X’,Y’).即:
$$X=OD,Y=AD$$
（2.1）求X’

从图中可以看出:
$$X^{\prime }=OE+BE$$
又因为
$$CD=BE$$
所以
$$X^{\prime }=OE+CD,OE=X\times \cos (\theta )$$
又因为
$$BE=CD=Y\times \sin (\theta )$$
最终得出：
$$X^{\prime }=OE+BE=X\times \cos (\theta )+Y\times \sin (\theta )\text{\hspace{1em}(1)}$$
（2.2）求Y’

由图可以很容易看出
$$Y^{\prime }=AC-BC$$

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qHlPLvEO-1650360946953)(C:\Users\pc\AppData\Roaming\Typora\typora-user-images\image-20220419144542247.png)]

又因为
$$AC=AD\times \cos (\theta ),BC=DE=OD\times \sin (\theta )$$
有上述所得:
$$AC=Y\times \cos (\theta ),BC=DE=X\times \sin (\theta )$$
最终得出：
$$Y^{\prime }=AB=AC-BC=Y\times \cos (\theta )-X\times \sin (\theta )\text{\hspace{1em}(2)}$$
由以上可以得出：
$$\{\begin{array}{l}{X}^{\prime }=X\times \cos (\theta )+Y\times \sin (\theta )\\ Y^{\prime }=Y\times \cos (\theta )-X\times \sin (\theta )\end{array}$$
对等式①与等式②同时提取
$$\cos (\theta )$$
得到：
$$\{\begin{array}{l}{X}^{\prime }=\cos (\theta )\times (X\times +Y\times \sin (\theta ))\\ Y^{\prime }=\cos (\theta )\times (Y\times -X\times \sin (\theta ))\end{array}$$
忽略cos()对于X’与Y’的影响，得到“伪旋公式”:
$$\{\begin{array}{l}{X}^{\prime }=X\times +Y\times \sin (\theta )\\ Y^{\prime }=Y-X\times \sin (\theta )\end{array}$$
伪旋公式：仅仅实现了正确的角度旋转，但向量模值变为了原来数值的
$$\frac{1}{\cos (\theta )}$$
使用CORDIC方法：

使得
$$\tan (\theta )=2^i，i\in -N^{\ast }$$

则“伪旋公式”进一步简化为：
$$\{\begin{array}{l}{X}^{\prime }=X\times +Y\times 2^i\\ Y^{\prime }=Y-X\times 2^i\end{array}$$
该式子已经仅剩下加减和移位运算了。

以下是
$$\tan (\theta )=2^i，i\in -N^{\ast }$$
所得的相关角度和正切值的对应表。

任意旋转角度θ，都可以由上表中的大小不同的θ进行多次旋转得到。CORDIC正是利用这一点，将θ角度的旋转分解为从大到小、逐次逼近真实旋转角度的一组旋转，而这些旋转的实现又都可以由加减和移位运算来完成。θi可以预先制成表格，供计算机查找使用。

（2）程序代码

#include <iostream>
#include <math.h>
#include <iomanip>　　　//设置必备的头文件
using namespace std;

void my_atan1(double A); //用于0~180 sin求解 
void my_atan2(double A); //用于180~360 sin求解
 
double S_C;

//主函数入口 
int main()
{	
    int	x;
    cout<<"enter a number:"<<endl;
	cin>>x;	
	//角度大于0 
	if(x>=0){
		x=x%360; 
	}
	//角度小于0 
	while(x<0){
		x+=360;
	}
	
	if(x>=0 && x<=90){
		my_atan1(x); //0~90
	}
	else if(x>90 && x<=180){
		my_atan1(180-x); //90~180
	}
	else if(x>180 && x<=270){
	    my_atan2(x-180);
	}else{
		my_atan2(360-x);
	}
	cout<<setiosflags(ios::fixed)<<setprecision(2);
	//cout<<"cos:"<< S_C[0] <<endl;
	cout<<"sin:"<<S_C<< endl;
	system("pause");
	return 0;
}

//角度属于0~180之间的角度数值求解 
void my_atan1(double A)
{	
	//迭代设置 
	double angle[]={45, 26.565, 14.036, 7.125, 3.5763,
                    1.7899, 0.8952, 0.4476, 0.22381,	
                    0.1119, 0.0559, 0.028, 0.01399		// 12 time
                   };
	double x=1;
	double y=0;
	double k=0.60723;//0.63664;
	int i=0;
	double x_new,y_new;
	//double angle=45;
	//进行迭代 
	for(i=0;i<12;i++)
	{
		angle[i]=(angle[i]/45)*pow(2,20);
	}
	A=(A/45)*pow(2,20);

	for(i=0;i<12; i++)
	{
		if(A>0)
		{
			x_new=(x-y*(1/pow(2,i)));
			y_new=(y+x*(1/pow(2,i)));
			x=x_new;
			y=y_new;
			A-=angle[i];
		}
		else{
			x_new=(x+y*(1/pow(2,i)));
			y_new=(y-x*(1/pow(2,i)));
			x=x_new;
			y=y_new;
			A+=angle[i];
		}
		//按对分查找方法寻找角度的近似值
		//angle/=2;
	}

	//S_C[0] = x * k;
	S_C=y*k;
}

//角度属于180~360之间的角度数值求解 
void my_atan2(double A)
{
	//迭代参数设置 
	double angle[]={45, 26.565, 14.036, 7.125, 3.5763,
                    1.7899, 0.8952, 0.4476, 0.22381,	
                    0.1119, 0.0559, 0.028, 0.01399		// 12 time
                   };
	double x=1;
	double y=0;
	double k=0.60723;//0.63664 0.60723;
	int i=0;
	double x_new,y_new;
	//double angle=45;
	//进行迭代 
	for(i=0;i<12;i++)
	{
		angle[i]=(angle[i]/45)*pow(2, 20);
	}
	A=(A/45)*pow(2,20);

	for(i=0;i<12;i++)
	{
		if(A>0)
		{
			x_new=(x-y*(1/pow(2,i)));
			y_new=(y+x*(1/pow(2,i)));
			x=x_new;
			y=y_new;
			A-=angle[i];
		}
		else{
			x_new=(x+y*(1/pow(2,i)));
			y_new=(y-x*(1/pow(2,i)));
			x=x_new;
			y=y_new;
			A+=angle[i];
		}
		//按对分查找方法寻找角度的近似值
		//angle/= 2;
	}

	//S_C[0]=x*k;
	S_C=-y*k;
}