原文转载自：https://zhuanlan.zhihu.com/p/20085048
在用C语言写程序的时候，如果需要某个数学方面的函数，包含一个math.h的头文件就万事大吉了。那么问题来了：如果没有这个math.h，用一些基本的数学知识，我们可以自己实现多少函数？
先看看math.h库里有哪些函数（注释是我加的）：
_CRTIMP double __cdecl sin (double);//正弦函数
_CRTIMP double __cdecl cos (double);//余弦函数
_CRTIMP double __cdecl tan (double);//正切函数
_CRTIMP double __cdecl sinh (double);//双曲正弦函数
_CRTIMP double __cdecl cosh (double);//双曲余弦函数
_CRTIMP double __cdecl tanh (double);//双曲正切函数
_CRTIMP double __cdecl asin (double);//反正弦函数
_CRTIMP double __cdecl acos (double);//反余弦函数
_CRTIMP double __cdecl atan (double);//反正切函数
_CRTIMP double __cdecl atan2 (double, double);//依然是反正切函数
_CRTIMP double __cdecl exp (double);//以e为底的指数函数
_CRTIMP double __cdecl log (double);//以e为底的对数函数
_CRTIMP double __cdecl log10 (double);//以10为底的对数函数
_CRTIMP double __cdecl pow (double, double);//实数的实数次方
_CRTIMP double __cdecl sqrt (double);//开平方
_CRTIMP double __cdecl ceil (double);//向正无穷方向取整
_CRTIMP double __cdecl floor (double);//向负无穷方向取整
_CRTIMP double __cdecl fabs (double);//绝对值
_CRTIMP double __cdecl ldexp (double, int);//计算前边那个实数乘e的后边那个整数次方的结果，没什么卵用
_CRTIMP double __cdecl frexp (double, int*);//把一个实数分解为x*2^exp次方的形式，貌似也没什么卵用
_CRTIMP double __cdecl modf (double, double*);//得到x的整数部分和小数部分
_CRTIMP double __cdecl fmod (double, double);//求实数除以实数的余数一般情况下，我们自己实现的函数在速度上肯定是不如库函数的，精度也很有限。因此本文选取的实现方式会适当地牺牲速度，但希望每个函数都可以达到1e-8的精度。下面我们就开始试着实现一些比较常用的函数。
1.基本常数
在函数实现的过程中，必不可少的要使用一些基本常数。我们首先对一些常用的函数进行宏定义，方便后边的程序编写。这些常数包括
[latex]\pi、\ln2、\ln10[/latex]
#define PI 3.14159265358979323846
#define e 2.7182818284590452354
#define ln_2 0.69314718055994530942
#define ln_10 2.302585092994045684022.取绝对值
为了更通用，我用宏定义来写绝对值函数：
#define fabs(a) ((a)>0?(a):(-(a)))3.实数的整数次方和实数的实数次方
首先，前者是非常好写的，直接用快速幂就可以。关于快速幂的介绍汗牛充栋，这里不再赘述。
double pow(double a,int n)
{
if(n<0) return 1/pow(a,-n);
double res = 1.0;
while(n)
{
if(n&1) res *= a;
a *= a;
n >>= 1;
}
return res;
}而实数的实数次方则需要使用一些数学知识了。由于[latex]a^{x}=\exp (x\ln a)[/latex],所以只要我们写出一个计算以[latex]e[/latex]为底的指数函数和以[latex]e[/latex]为底的对数函数，该函数自然就出来了。
double powf(double a,double x)
{
return exp(x*ln(a));
}在这里留下了两个坑：[latex]\exp(x)、\ln x[/latex]下边将分别填坑。
4.以e为底的指数函数
考虑到[latex]f(x)=\exp(x)[/latex]定义域为整个实数域，因此有必要对区间进行一些调整。我们先令[latex]x=0[/latex]时将[latex]x[/latex]取相反数，同时整个函数值取倒数，即
if(x < 0) return 1 / exp(-x);然后将[latex]x[/latex]分为整数部分和小数部分，整数部分用实数的整数次方的函数计算，小数部分单独计算。代码如下：
int n = (int)x;
x -= n;
double e1 = pow(e,n);
double e2 = eee(x);
return e1*e2;那个eee就是计算小数部分的函数。我们试着用泰勒展开式直接计算：
double eee(double x)
{
return 1 + x + x*x/2 + pow(x,3)/6 + pow(x,4)/24 + pow(x,5)/120;
}用卡西欧991计算器和该程序分别计算[latex]\exp(5.3333)[/latex]，得数分别为207.1203448和207.120048，精度显然不够。怎么办呢？因为该展开式是在[latex]x=0[/latex]这一点的，越靠近0就越精确。而[latex] \exp(x)=\exp^{2}(\frac{x}{2}) [/latex],故我们可以递归地进行这样的代换，直到[latex]x[/latex]足够小。即：
double eee(double x)
{
if(x>1e-3)
{
double ee = eee(x/2);
return ee*ee;
}
return 1 + x + x*x/2 + pow(x,3)/6 + pow(x,4)/24 + pow(x,5)/120;
}这时计算结果为207.12034476，精度足够了！
5.自然对数和任意底数对数
因为
[latex]\Large \ln x=\int_{1}^{x}\frac{dt}{t}[/latex]
所以用自适应辛普森公式计算这个积分即可。辛普森公式：
[latex]\Large \int_{a}^{b} f(x)dx \approx \frac{b-a}{6} \left[ f(a)+4f\left( \frac{a+b}{2} \right) +f(b) \right][/latex]
double F1(double x)
{
return 1/x;
}
double simpson(double a, double b,int flag)
{
double c = a + (b-a)/2;
if(flag==1)
return (F1(a)+4*F1(c)+F1(b))*(b-a)/6;
if(flag==2)
return (F2(a)+4*F2(c)+F2(b))*(b-a)/6;
}
double asr(double a, double b, double eps, double A,int flag)
{
double c = a + (b-a)/2;
double L = simpson(a, c,flag), R = simpson(c, b,flag);
if(fabs(L+R-A) <= 15*eps) return L+R+(L+R-A)/15.0;
return asr(a, c, eps/2, L,flag) + asr(c, b, eps/2, R,flag);
}
double asr(double a, double b, double eps,int flag)
{
return asr(a, b, eps, simpson(a, b,flag),flag);
}
double ln(double x)
{
return asr(1,x,1e-8,1);
}而任意底数的对数非常简单，应用换底公式[latex]\log_{a} x=\frac{\ln b}{\ln a}[/latex]即可。
double log(double a,double N)
{
return ln(N)/ln(a);
}6.开平方算法
大家一定都知道卡马克那个超级快的平方根倒数算法，但是我不打算那样做，而是用最土鳖的牛顿迭代：给定初值[latex]a_{0}=x/2[/latex]（或其他）,然后用递推式
[latex]\Large a_{n+1}=\frac{1}{2} \left( a_{n}+\frac{x}{a_{n}} \right)[/latex]
计算。代码如下（我先迭代了两次，把迭代结果作为初值）：
double sqrt(double x)
{
double t = x/8 + 0.5 + 2*x/(4+x);
int c = 10;
while(c--)
{
t = (t+x/t)/2;
}
return t;
}这里我迭代十次。精度如何呢？随便试几个数：
[latex]\sqrt{98756} = 314.25467379[/latex]，和计算器一模一样，完美。
[latex]\sqrt{135791113} = 19220.7977[/latex]，然而真实值为11652.944，差得太远了！
当然也可以用自适应精度，但是联想到计算[latex][/latex]为底的指数函数时采用的方法，我们可以令一个数大于100时不停地除以100，同时给结果乘10.这样最后的数就是100以内的，精度会很高。而一个常见的数除以100是除不了几次的。代码如下：
double sqrt(double x)
{
if(x>100) return 10.0*sqrt(x/100);
double t = x/8 + 0.5 + 2*x/(4+x);
int c = 10;
while(c--)
{
t = (t+x/t)/2;
}
return t;
}这下结果完美了，再大几千倍的数也不成问题。
7.三角函数
三角函数的计算我在这个回答里介绍了，同样是缩小区间再泰勒。
利用周期性把[latex]x[/latex]变成[latex](-\pi ,\pi)[/latex]中的数，
再利用奇偶性把[latex]x[/latex]变成[latex][0,\pi][/latex]中的数，
再用诱导公式把[latex]x[/latex]变成[latex][0,\frac{\pi}{2}][/latex]中的数，
用诱导公式把[latex]x[/latex]变成[latex][0,\frac{\pi}{4}][/latex]中的数
最后级数展开。
double sin(double x)
{
double fl = 1;
if(x>2*PI || x<-2*PI) x -= (int)(x/(2*PI))*2*PI;
if(x>PI) x -= 2*PI;
if(x<-PI) x += 2*PI;
if(x>PI/2)
{
x -= PI;
fl *= -1;
}
if(x<-PI/2)
{
x += PI;
fl *= -1;
}
if(x>PI/4) return cos(PI/2-x);
else return fl*(x - pow(x,3)/6 + pow(x,5)/120 - pow(x,7)/5040 +pow(x,9)/362880);
}
double cos(double x)
{
double fl = 1;
if(x>2*PI || x<-2*PI) x -= (int)(x/(2*PI))*2*PI;
if(x>PI) x -= 2*PI;
if(x<-PI) x += 2*PI;
if(x>PI/2)
{
x -= PI;
fl *= -1;
}
if(x<-PI/2)
{
x += PI;
fl *= -1;
}
if(x>PI/4) return sin(PI/2-x);
else return fl*(1 - pow(x,2)/2 + pow(x,4)/24 - pow(x,6)/720 + pow(x,8)/40320);
}正切函数最简单：两个一除完事。
double tan(double x)
{
return sin(x)/cos(x);
}这种方法的精度居然也是足够的。
8.反三角函数
反三角函数用级数计算的效果似乎不太好，不知道为什么。然而出乎意料的是，用定积分精度却很不错。当然，这个精度依赖于我们写的开平方函数的精度。原理为
[latex]\Large \arcsin (x)=\int_{0}^{x} \frac{dt}{\sqrt{1-t^{2}}}[/latex]
代码如下（特判一下正负1，因为这两点分母为0）：
double F2(double x)
{
return 1/sqrt(1-x*x);
}
double simpson(double a, double b,int flag)
{
double c = a + (b-a)/2;
if(flag==1)
return (F1(a)+4*F1(c)+F1(b))*(b-a)/6;
if(flag==2)
return (F2(a)+4*F2(c)+F2(b))*(b-a)/6;
}
double asr(double a, double b, double eps, double A,int flag)
{
double c = a + (b-a)/2;
double L = simpson(a, c,flag), R = simpson(c, b,flag);
if(fabs(L+R-A) <= 15*eps) return L+R+(L+R-A)/15.0;
return asr(a, c, eps/2, L,flag) + asr(c, b, eps/2, R,flag);
}
double asr(double a, double b, double eps,int flag)
{
return asr(a, b, eps, simpson(a, b,flag),flag);
}
double asin(double x)
{
if(fabs(x)>1) return -1;
double fl = 1.0;
if(x<0) {fl*=-1;x*=-1;}
if(fabs(x-1)<1e-7) return PI/2;
return (fl*asr(0,x,1e-8,2));
//return x + pow(x,3)/6 + pow(x,5)*3/40 +pow(x,7)*5/112 + pow(x,9)*35/1152 + pow(x,11)*315/1408;
}经试验，精度也是完美的。
反余弦函数直接用[latex]\frac{\pi}{2}[/latex]减就行：
double acos(double x)
{
if(fabs(x)>1) return -1;
return PI/2 - asin(x);
}而反正切函数，直接用泰勒展开到9次方精度依然不足，区间调整到[latex][0,1][/latex]之后精度有所提高， 但依然只有5位左右。看来只有用一点数学技巧了。
若[latex]xy<1[/latex]，有
[latex]\large \arctan x+\arctan y=\arctan \frac{x+y}{1-xy}[/latex]
当[latex]x=y[/latex]时，该式就变成了
[latex]\large 2\arctan x=\arctan \frac{2x}{1-x^{2}}[/latex]
令[latex]\large \frac{2x}{1-x^{2}} =y[/latex]，解得[latex]\large x=\frac{-1+\sqrt{1+y^{2}}}{y}[/latex]
(不合理根已舍去)
这个值是小于[latex]y[/latex]的，故可以用这个式子迭代使自变量足够小以满足泰勒公式的精度。
double atan(double x)
{
if(x<0) return -atan(-x);
if(x>1) return PI/2 - atan(1/x);
if(x>1e-3) return 2*atan((sqrt(1+x*x)-1)/x);
return x - pow(x,3)/3 + pow(x,5)/5 - pow(x,7)/7 + pow(x,9)/9;
}9.双曲函数
因为以[latex]e[/latex]为底的指数函数已经有了，所以双曲函数就不用啰嗦了。直接上代码：
double sinh(double x)
{
return (exp(x)-exp(-x))/2;
}
double cosh(double x)
{
return (exp(x) + exp(-x))/2;
}
double tanh(double x)
{
double ext = exp(x);
return 1 - 2/(ext*ext+1);
}10.总结
至此，我们自己实现了math.h中的大多数函数，速度可能逊色一些，但精度够用。最关键的是，这些实现方式采用的数学方法都是很简单易懂的，理解之后完全可以自己写出来。
11.最终代码（加注释）
#include <stdio.h>
#define PI 3.14159265358979323846
#define e 2.7182818284590452354
#define ln_2 0.69314718055994530942
#define ln_10 2.30258509299404568402
#define fabs(a) ((a)>0?(a):(-(a)))
#define maxf(a,b) ((a) > (b) ? (a) : (b))
#define minf(a,b) ((a) < (b) ? (a) : (b))
double pow(double a,int n);
double powf(double a,double x);
double sqrt(double x);
double exp(double x);
double ln(double x);
double cos(double x);
double sin(double x);
double tan(double x);
double simpson(double a, double b,int flag);
double asr(double a, double b, double eps, double A,int flag);
double asr(double a, double b, double eps,int flag);
double acos(double x);
double asin(double x);
double atan(double x);
double cosh(double x);
double sinh(double x);
double tanh(double x);
double log(double a,double N);//函数声明
double pow(double a,int n)//实数的整数次幂，快速幂
{
if(n<0) return 1/pow(a,-n);//n为负时的处理
double res = 1.0;
while(n)
{
if(n&1) res *= a;
a *= a;
n >>= 1;
}
return res;
}
double powf(double a,double x)//实数的实数次幂，利用exp函数
{
return exp(x*ln(a));
}
double F1(double x)//为了使用定积分计算lnx
{
return 1/x;
}
double F2(double x)//为了使用定积分计算arcsinx
{
return 1/sqrt(1-x*x);
}
double simpson(double a, double b,int flag)//辛普森积分公式
{
double c = a + (b-a)/2;
if(flag==1)
return (F1(a)+4*F1(c)+F1(b))*(b-a)/6;
if(flag==2)
return (F2(a)+4*F2(c)+F2(b))*(b-a)/6;
}
double asr(double a, double b, double eps, double A,int flag)//辛普森积分公式的自适应过程
{
double c = a + (b-a)/2;
double L = simpson(a, c,flag), R = simpson(c, b,flag);
if(fabs(L+R-A) <= 15*eps) return L+R+(L+R-A)/15.0;
return asr(a, c, eps/2, L,flag) + asr(c, b, eps/2, R,flag);
}
double asr(double a, double b, double eps,int flag)//辛普森积分主体函数
{
return asr(a, b, eps, simpson(a, b,flag),flag);
}
double sin(double x)
{
double fl = 1;
if(x>2*PI || x<-2*PI) x -= (int)(x/(2*PI))*2*PI;
if(x>PI) x -= 2*PI;
if(x<-PI) x += 2*PI;
if(x>PI/2)
{
x -= PI;
fl *= -1;
}
if(x<-PI/2)
{
x += PI;
fl *= -1;
}
if(x<0)
{
x*=-1;
fl*=-1;
}//处理区间
if(x>PI/4) return cos(PI/2-x);
else return fl*(x - pow(x,3)/6 + pow(x,5)/120 - pow(x,7)/5040 +pow(x,9)/362880);//泰勒公式
}
double cos(double x)
{
double fl = 1;
if(x>2*PI || x<-2*PI) x -= (int)(x/(2*PI))*2*PI;
if(x>PI) x -= 2*PI;
if(x<-PI) x += 2*PI;
if(x>PI/2)
{
x -= PI;
fl *= -1;
}
if(x<-PI/2)
{
x += PI;
fl *= -1;
}//处理区间
if(x>PI/4) return sin(PI/2-x);
else return fl*(1 - pow(x,2)/2 + pow(x,4)/24 - pow(x,6)/720 + pow(x,8)/40320);//泰勒公式
}
double tan(double x)
{
if(x>PI || x<-PI) x -= (int)(x/(PI))*PI;
if(x>PI/2) x -= PI;
if(x<-PI/2) x += PI;//处理区间
return x + pow(x,3)/3 + pow(x,5)*2/15 + pow(x,7)*17/315 + pow(x,9)*62/2835;//泰勒公式
}
double asin(double x)
{
if(fabs(x)>1) return -1;
double fl = 1.0;
if(x<0) {fl*=-1;x*=-1;}
if(fabs(x-1)<1e-7) return PI/2;//x为正负1时特判
return (fl*asr(0,x,1e-8,2));//主体积分过程
//return x + pow(x,3)/6 + pow(x,5)*3/40 +pow(x,7)*5/112 + pow(x,9)*35/1152 + pow(x,11)*315/1408;
}
double acos(double x)
{
if(fabs(x)>1) return -1;
return PI/2 - asin(x);//简单公式的应用
}
double atan(double x)
{
if(x<0) return -atan(-x);
if(x>1) return PI/2 - atan(1/x);
if(x>1e-3) return 2*atan((sqrt(1+x*x)-1)/x);//递推地缩小自变量，使之接近0，保证泰勒公式的精度
return x - pow(x,3)/3 + pow(x,5)/5 - pow(x,7)/7 + pow(x,9)/9;//泰勒公式
}
double eee(double x)//计算exp函数中的小数部分
{
if(x>1e-3)
{
double ee = eee(x/2);//递推缩小自变量，使之接近0，保证精度
return ee*ee;
}
else
return 1 + x + x*x/2 + pow(x,3)/6 + pow(x,4)/24 + pow(x,5)/120;
}
double exp(double x)
{
if(x<0) return 1/exp(-x);
int n = (int)x;
x -= n;//分割整数部分和小数部分分别计算
double e1 = pow(e,n);//计算整数部分
double e2 = eee(x);//计算小数部分
return e1*e2;//计算结果
}
double sinh(double x)
{
return (exp(x)-exp(-x))/2;
}
double cosh(double x)
{
return (exp(x) + exp(-x))/2;
}
double tanh(double x)
{
double ext = exp(x);
return 1 - 2/(ext*ext+1);
}
double ln(double x)
{
return asr(1,x,1e-8,1);//直接算积分
}
double sqrt(double x)
{
if(x>100) return 10.0*sqrt(x/100);//缩小自变量值
double t = x/8 + 0.5 + 2*x/(4+x);//迭代两次之后的结果
int c = 10;//迭代十次
while(c--)
{
t = (t+x/t)/2;
}
return t;
}
double hypot(double x,double y)
{
return sqrt(x*x + y*y);//这是一个已知直角三角形两直角边求斜边的程序，没有什么卵用
}
double log(double a,double N)
{
return ln(N)/ln(a);//换底公式计算任意对数
}
int main()
{
printf("%.10lf\n",atan(0.3));
return 0;
}