//////////////////////////////////////////////////////////////////////////
//Calculat logarithm
//Aitken算法加速數列收斂算法
double expP(double x)//計算e^x
{
double y = x;
double ex_p1 = 0;
double ex_p2 = 0;
double ex_p3 = 0;
double ex_p = 0;
double ex_px = 0;
double ex_tmp = 1;
double dex_px = 1;
double tmp;
int l;
if (0 == x)
{
return 1;
}
if(x < 0)
{
return 1/expP(-x);
}
for(l = 1,tmp = 1;(ex_px - ex_tmp)>0.0000000001 || (ex_px-ex_tmp)<(-0.0000000001) && dex_px > 0.0000000001;l++)
{
ex_tmp = ex_px;
tmp *= y;
tmp = tmp/l;
ex_p1 += tmp;
ex_p2 = ex_p1 + tmp * y / (l+1);
ex_p3 = ex_p2 + tmp * y * y / (l+1) /(l+2);
dex_px = ex_p3 - ex_p2;
ex_px = ex_p3 - dex_px *dex_px / (ex_p3 - 2 * ex_p2 + ex_p1);
}
return ex_px + 1;
}
//////////////////////////////////////////////////////////////////////////
double InP(double x)//計算In(x)
{
double y = x -1;
double In_p1 = 0;
double In_p2 = 0;
double In_p3 = 0;
double In_p = 0;
double In_px = 0;
double In_tmp = 1;
double dIn_px = 1;
double tmp;
int l;
if(1 == x)
{
return 0;
}
else if(x > 2)
{
tmp = -InP(1/x);
return tmp;
}
else if(x < 1)
{
double n = -1;
double a;
do
{
n = n - 0.06;
a = x / expP(n);
}
while (a > 2|| a < 1);
tmp = InP(a) + n;
return tmp;
}
for(l = 1, tmp = 1.00;(In_px - In_tmp) > 0.0000000001 || (In_px - In_tmp) < -0.0000000001; l++)
{
In_tmp = In_px;
tmp *= y;
if(1== l)
{
tmp = tmp / l;
}
else
{
tmp = tmp / (-l);
}
In_p1 += tmp;
In_p2 = In_p1 +(-1) * tmp * y * l / (l + 1);
In_p3 = In_p2 + tmp * y * y * l / (l + 2);
dIn_px = In_p3 -In_p2;
In_px = In_p3 - dIn_px * dIn_px / (In_p3 - 2 * In_p2 + In_p1);
tmp *= l;
}
return In_px;
}
1 unsigned int fn(int n)
2 {
3 return (n == 0 || n == 1)? 1 : n * fn(n-1);
4 }
5
6 double Mysin(double x)
7 {
8 int m = 1, sign = 1;
9 double t, sum = 0;
10
11 while ( fabs(t = sign * pow(x, 2*m -1) / fn(2*m - 1)) > 0.00001)
12 {
13 sum += t;
14 ++m;
15 sign *= -1;
16 printf("t=%lf sum = %lf\r\n",t,sum);
17 }
18
19 return sum;
20 }
21 //////////////////////////////////////////////////////////////////////////
22 double Mycos(double y)
23 {
24 return (sqrt(1-Mysin(y)*Mysin(y)));
25 }
26 //////////////////////////////////////////////////////////////////////////
27