C++代码实现atanx,atan2x,acosx,asinx,sinx,cosx,tanx三角函数、开方函数sqrtND()

一、开平方函数，牛顿迭代法sqrtND()

double sqrtND(double A) //牛顿迭代法
{
double k = A;
while (kk - A > 0.00000001 || k * k - A < -0.00000001)
{
k = (0.5k + A / 2 / k) ;
}
return k;
}

二、sinx()/cosx()/tanx()

double angle[] = { pi / 4,pi / 8,pi / 16,
pi / 32,pi / 64,pi / 128,
pi / 256,pi / 512,pi / 1024,
pi / 2048,pi / 4096,pi / 8192,pi / 16384 };

double tang[] = { 1,0.4142135623731,0.19891236737966,
0.098491403357164,0.049126849769467,
0.024548622108925,0.012272462379566,
0.0061360001576234,0.0030679712014227,
0.0015339819910887,0.00076699054434309,
0.00038349521577144,0.00019174760083571 };

double sinx(double a) {
while (a < 0)a += 2 * pi;
if (a <= (pi / 16384)) {
return a; //因为a的值太小，sin a 约等于 a
}
else {
//开始 CORDIC 算法
bool negitive = a < 0; // sin(-a) = -sin(a)
double x = 10;
double y = 0;
double theta = 0;
for (int i = 0;i < 13;i++) { //开始迭代
double orix=0, oriy=0;
while (theta < a) { //当前角度小于a
orix = x;
oriy = y;
//坐标旋转
x = orix - tang[i] * oriy;
y = tang[i] * orix + oriy;
theta += angle[i];
}
if (theta == a) {
if (negitive) {
return -(y / sqrtND((xx + y * y)));
}
else {
return (y / sqrtND((xx + y * y)));
}
}
else {
//旋转的弧度超过了a，退回原来增加的角度，同时进入下一次迭代
theta -= angle[i];
x = orix;
y = oriy;
}
}
if (negitive) {
return -(y / sqrtND((xx + y * y)));
}
else {
return (y / sqrtND((xx + y * y)));
}
}
}

double cosx(double a)
{
return sinx(pi/2-a);
}

三、atanx()/asinx()/acosx()/atan2x()

/三角函数全部用泰勒展开，其中atanx收敛域为-1~1，需进行转化/
double atanx(double a)
{
/值域[-pi/2,pi/2]/
/泰勒展开在-1~1内是收敛的，超出范围发散，可进行转换1/x，x>1,atan(x)=pi/2-atan(1/x),x<-1,atan(x)=-pi/2-atan(1/x)/
if (a == 1)return pi / 4;
else if (a == -1)return -pi / 4;
if (a<1 && a>-1)
{
double p = 0.000000001, g = 0, t = a, tmp;
int n = 0;
do {
g = g + t / (2 * n + 1);
t = -t * a*a;
n++;
tmp = t / (2 * n + 1);
} while (fabs(tmp) > p);
return g;
}
else if(a<-1)
return -pi_2-atanx(1 / a);
else return pi_2 - atanx(1 / a);

}
double asinx(double a)
{
/值域[-pi/2,pi/2]/
if (a == 0.5)return pi / 6;
else if (a == -0.5)return -pi / 6;

double p = 0.00000001, g = a, t = a, tmp=a;
int n = 1;
do {tmp = tmp* (2*n-1)*(2*n-1) * a*a / (2*n) / (2*n + 1);g = g + tmp;
} while (fabs(tmp) > p);
return g;

}
double acosx(double a)
{
/值域[0,pi]/
return pi_2 - asinx(a);
}
double atan2x(double y, double x)
{
/值域[-pi,pi]/
if (x > 0) return atanx(y / x);
else if (x < 0 && y >= 0)
return atanx(y / x) + pi;
else if (x < 0 && y < 0)
return atanx(y / x) - pi;
else if (x == 0 && y > 0)
return pi_2;
else if (x == 0 && y < 0)
return -pi_2;
else return 0.0;
}