最近接触到CORDIC算法，说是一种使用迭代去计算三角函数的做法。我看了原理，感觉写的晦涩难懂，但本质上是先计算一个初始角度，然后再慢慢增加角度，从而得出大角度的值。但是我转念一想，既然大角度的值可以分解为小角度来计算，我何必搞得这么复杂呢？
我们知道
sin(A+B) = sin(A)*cos(B)+cos(A)*sin(B);
cos(A+B) = cos(A)*cos(B)-sin(A)*sin(B);
假设我需要计算sin(30).那么由于:
sin(2)=sin(1+1)= sin(1)*cos(1)+cos(1)*sin(1);
sin(4)=sin(2+2)= sin(2)*cos(2)+cos(2)*sin(2);
.....
sin(30)=sin(15 + 15) =sin(15)*cos(15)+cos(15)*sin(15);
所以，我可以先计算出sin(1),cos(1)就能得到sin(2),cos(2);然后根据sin(2),cos(2)又能计算出sin(4),cos(4)。。。最终我能计算出sin(30)。换句话说，只要我知道sin(1),cos(1)的值，那么0~90度的sin值我都能算。而360度都能转化到0~90度之间来算。当然,sin(0)和cos(0)的值是显而易见的，无须计算。
那么:只要我知道sin(1),cos(1)的值，我就能逐步推算出0~360度任意一个整数角度的三角函数值。
而sin(1) = 0.01745240644,cos(1) = 0.99984769516，可以直接定义成常量。
一个显著的缺点是，角度的分辨率是1度。如果要提升，可以把原始值sin(1),cos(1)换成更小的角度即可。比如计算sin(30)，角度分割的过程如下图所示，最终分解成只计算1,0的三角函数。
这里总结一下，在没有标准数学库的情况下，自己实现三角函数计算的常用方法。
一，查表法
直接把0~90的角度按固定的步长取值，把对应的sin事先算好，存储在一个数组里面。在需要计算sin(x)时，直接取与x最接近的点查表即可。
这种方法速度快，编程简单，但精度一高就要占用很多存储空间。在不需要计算精度多高的情况下，不失为一个好选择。
二，尊敬的小蝴蝶级数展开
根据尊敬的小蝴蝶公式，sinx可以展开为无穷级数sinx=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!-……。按这个公式进行累加运算。
这种方法需要计算很多乘法、除法、计算量较大。但随着迭代次数增多，精度越高。且精度可控，当达到需要的精度时可停止计算。但由于计算复杂，可能输出结果的速度较慢。
三，递归分解法
就像本文提到的，把大角度一分为二，再分，直至为0或1来计算。
这种方法可以达到很高精度，编程也简单，但是不断的分割角度需要函数递归。当设置的精度越高，函数递归的层数显著增多。需要注意栈空间不足的问题。
在计算大角度时，由于递归层数太多，容易遇到计算速度太慢的问题。
同时也需要乘法运算。
四，cordic算法
感觉cordic算法和第三种方法根本就是一回事.假设A为原始角度，B为旋转角度，那么旋转之后的sin值。
sin(A+B) = sin(A)*cos(B)+cos(A)*sin(B);
然后取B=45,30,15,5,2,1这样的特殊值，事先存储在代码中，就能快速计算。
比如要算sin(75)，可设置原始角度A=0.
然后计算sin(0+50)=>sin(50)=>sin(50+30)=>sin(80)=>sin(80-5)=>sin(75)。
而cordic算法提到的"只需要使用移位和加减法"，根本就是个噱头。因为乘法本来就可以转换成移位和加减法，和cordic毫无关系。它只是选择了一些特殊的角度，让计算出来的值刚好适合移位运算而已。比如sin(0.01398822727)=2^-12,它就选这种0.01398822727为旋转角度。当数字X乘以sin(0.01398822727)时，就右移动12bit即可。
当然，尽管把cordic说的这么简单，它的优点还是要学习的，就是不像第三种方法一样，只能分解角度到1才能计算。我们可以学习这一点，来改善方法三的代码。即:事先将sin(60),sin(45),sin(10),sin(5),sin(2),sin(1)等特殊值存储在数组中，角度分解时可以尽量往特殊角度上靠，实在不能分解的才分解为1。为了可以计算小数角度，规定在表中的角度为扩大100倍之后的值。这样，三角函数计算角度可以精确到小数点后两位。修改后的示意图如下。可见比起方法3，程序执行递归的次数远远减少了。比如计算sin(77)，角度分割过程如下:
要注意这个sin_table表。这个表里面存储的值越多，函数递归的次数原则上是会减少的，也就是程序运行的更快。但存的多占用的空间也大，所以一般选90度，45度，30度，10度，5度，2度，1度，0度的值。要注意的是，这个表就算再小，至少也要存储0度的值和1度的值。如果要求角度精确到0.1度，那要存储0.1度的三角函数值。也就说0度和最小角度必须存储。这是因为角度最小只能分解为0和最小角度。
最小角度决定了计算角度的分辨率，表中存储数据多少影响了计算速度。
但并不是表中存储数据越多就一定能加快速度。如果是[92度，91度，90度，1度，0度]这样的表，起不到明显的加速作用。必须合理分配角度值。
但是注意，这并不是真正的cordic分解角度的流程。cordic分解角度，是每一步都是加上或者减去固定的角度来的。它不再考虑角度应该怎么分解成两个什么小角度，因为每一步变化的角度值是固定的，不同的只是选择加上或者减去。可见下表。
假设计算sin(41),初始角度0.角度逼近过程是:
第i步当前角度cordic角度表决策逼近结果0045+4514526.56505118-18.43494882218.4349488214.03624347+32.47119229332.471192297.125016349+39.59620864439.596208643.576334375+43.17254301543.172543011.789910608-41.38263241641.382632410.89517371-40.4874587740.48745870.447614171+40.93507287840.935072870.2238105+41.15888337941.158883370.111905677-41.046977691041.046977690.005595289-41.0413824
这样做的好处就是，由于每一步变化的角度是固定的，所以可以从sin(A+B) = sin(A)*cos(B)+cos(A)*sin(B)中又继续提取一个公因子cos(B)，来简化运算。由于cos函数是偶函数，所以cos(-B)=cos(B)。角度逐步逼近时，不论你每一步是加还是减，都能视为cos(B)。假设分解角度一共进行了i步，那么K=cos(B0) *cos(B1)*...*cos(Bi)也可以查表，所以无需每一步都计算cos(B),只需要记录步数，然后查表算出最终K即可。
但这样做的坏处就是，明明一个很容易分解的角度，非要加来减去的分解很多次。按之前的说法，41度只需分解为30度和11度，继而将11度分解为10度和1度即可。
结合了第三种方法和cordic优点的代码如下：
#include <stdio.h>
#define ANGLE_ZOOM 100
double my_int_sin(int x);
double my_int_cos(int x);
typedef struct SIN_ITEM
{
int angle;
double sin_value;
double cos_value;
}SIN_ITEM;
//角度值扩大100倍存储
SIN_ITEM sin_table[] =
{
9000, 1, 0,
4500, 0.70710678119, 0.70710678119,
3000, 0.5, 0.86602540378,
500, 0.08715574275, 0.99619469809,
200, 0.0348994967, 0.99939082702,
100, 0.01745240644, 0.99984769516,
50, 0.0087265355, 0.99996192306,
20, 0.00349065142, 0.99999390766,
10, 0.00174532837, 0.99999847691,
5, 0.00087266452, 0.99999961923,
2, 0.00034906584, 0.99999993908,
1, 0.00017453292, 0.99999998477,
0, 0, 1,
};
double my_int_sin(int x)
{
int i = 0;
int A=0,B=0;
for(i=0;i<sizeof(sin_table)/sizeof(sin_table[0]);i++)
{
if(x == sin_table[i].angle)
{
return sin_table[i].sin_value;
}
}
for(i=0;i<sizeof(sin_table)/sizeof(sin_table[0]);i++)
{
if(sin_table[i].angle <= x && sin_table[i-1].angle> x)
{
A= sin_table[i].angle;
}
}
if(A==0)
{
A= sin_table[0].angle;
}
B = x - A;
return my_int_sin(A)*my_int_cos(B)+my_int_cos(A)*my_int_sin(B);
}
double my_int_cos(int x)
{
int i = 0;
int A=0,B=0;
for(i=0;i<sizeof(sin_table)/sizeof(sin_table[0]);i++)
{
if(x == sin_table[i].angle)
{
return sin_table[i].cos_value;
}
}
for(i=0;i<sizeof(sin_table)/sizeof(sin_table[0]);i++)
{
if(sin_table[i].angle <= x && sin_table[i-1].angle> x)
{
A= sin_table[i].angle;
}
}
if(A==0)
{
A= sin_table[0].angle;
}
B = x - A;
return my_int_cos(A)*my_int_cos(B)-my_int_sin(A)*my_int_sin(B);
}
double my_sin(double x)
{
int angle = x*ANGLE_ZOOM;
int sign =1;
angle = angle % (360*ANGLE_ZOOM);
if(angle < 0)
{
sign = -1;
angle = -angle;
}
return sign*my_int_sin(angle);
}
double my_cos(double x)
{
int angle = x*ANGLE_ZOOM;
int sign =1;
angle = angle % (360*ANGLE_ZOOM);
if(angle < 0)
{
angle = -angle;
}
return sign*my_int_cos(angle);
}
void main()
{
double x = 0;
x = -43.97;
printf("angle %f sin= %.10f cos= %.10fn",x,my_sin(x),my_cos(x));
x = -598.32;
printf("angle %f sin= %.10f cos= %.10fn",x,my_sin(x),my_cos(x));
}
运行结果，可见基本正确: