数值算法(一)——自然常数为底的对数(ln)和指数函数(exp)

之前在使用SIMD的时候发现Intel官网上的指南中是有log、pow等函数的，但是用gcc和clang都不能编译，查了一下才知道Intel的CPU只支持四则运算和sqrt五个计算指令，而pow、log这类函数icc是通过库函数实现的。为了能用gcc编译不得已只能考虑自己写一个，本来以为很简单的任务，实际上还确实没那么好写。

exp函数的实现

对于exp(x)函数实现计算$e^x$ex，对$e^x$ex在$x=0$x=0处做泰勒展开可以得到
$e^x=1+x+\frac{1}{2}x^2+\frac{1}{3}x^3+.....+\frac{1}{n!}x^n$ex=1+x+21​x2+31​x3+.....+n!1​xn
利用这个多项式即可实现exp函数。这里需要考虑两个问题，第一这个级数的收敛条件、第二截断误差。

这个级数的收敛条件很简单，只要x是有限值，级数收敛。而截断误差就比较麻烦了，因为通常我们希望截断误差应该是一个固定的值，但是这个式子计算的截断误差是一个与x和n都有关的值。我们当然可以从第一项开始通过一个while循环来判断每一次计算与上一次结果的差值，从而决定是否继续计算下去。

如果从1开始一项一项向后加，加到这个多项式的第n项则需要$O(n^2)$O(n2)次乘法和$O(n)$O(n)次加法。如果使用秦九韶算法求这个多项式则只需要$O(n)$O(n)次乘法和加法。但是秦九韶算法需要我们从最高阶项的系数开始计算，而如果我们想保证截断误差固定，我们就需要每次计算前先计算n的值，然后从第n项的系数开始计算，这很繁琐。并且如果x很大n也会比较大，这时需要的计算量也会增加，并不高效。

于是我们考虑把x拆分为整数部分i和小数部分m两部分。
$x=i+m$x=i+m
对于整数部分i我们可以直接利用i个$e$e相乘直接得到，而且还可以使用快速幂算法计算。而小数部分m必然满足$0<m<1$0<m<1，这部分用泰勒展开近似，由于$m<1$m<1所以假设我们需要保证$10^{-14}$10−14的截断误差，我们只需要最高计算到16项即可，这个计算量相对来说并不算多，因此可以不考虑m的实际取值，而都计算到第16项截断，必然可以满足相应的截断误差。

为了进一步减小计算量，我们对m再进行一次处理。可以看到如果截断误差固定，那么m的值越小，需要的阶数就越低。那么我们直接用m除以2的幂次，来减小m的值，最后只需要把计算再不断平方就能算回原来的值（$e^{m}=(e^{m/2^{k}})^{2^{k}}$em=(em/2k)2k）。在计算$e^{m/2^{k}}$em/2k时我们相当于把m缩小了，于是可以使用更低的阶数达到更高的精度，这节约了计算量，而在计算$e^{m}=(e^{m/2^{k}})^{2^{k}}$em=(em/2k)2k时，通过快速幂算法，使得原本$O(2^{k})$O(2k)次的乘法变成了$O(k)$O(k)次乘法，这样取适当的k值可以显著减小总计算量。这里没有做过多的测试，直接取k为4，截断误差取$10^{-14}$10−14时，需要阶数为7。于是得到下面的代码。

#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;
const double _e=2.718281828459045;
double _exp_i(const int &x)
{
	double r=1,tmp=_e;
	int t=x;
	while(t)
	{
		if(t%2) r*=tmp;
		tmp*=tmp;
		t/=2;
	}
	return r;
}
double _exp_d(const double &x)
{
	
	double r=0,dx=x/16,f=1.0/5040;
	for(int i=7;i>0;i--)
	{
		r+=f;
		r*=dx;
		f*=i;
	}
	r+=1;
	r*=r;
	r*=r;
	r*=r;
	return r*r;
}
double _exp(const double &x)
{
	int i=round(x);
	double r=_exp_d(x-i);
	return r*_exp_i(i);
}
int main()
{
	double dt;
	double i=double(clock()%16)/16;
	cout<<"error = "<<exp(i)-_exp(i)<<endl;
	clock_t s=clock();
	double res[1000000];	
	for(int i=1;i<1000000;i++)
	{
		res[i]=exp(double(i)/50000);
	}
	dt=-double(s-clock())/CLOCKS_PER_SEC;
	cout<<res[7777]<<endl;
	s=clock();
	for(int i=1;i<1000000;i++)
	{
		res[i]=_exp(double(i)/50000);
	}
	cout<<"_exp "<<-double(s-clock())/CLOCKS_PER_SEC<<endl;
	cout<<"exp "<<dt<<endl;
	cout<<res[7777]<<endl;
	return 0;
}

使用g++ -O3编译运行得到

_exp 0.020947
exp 0.024846

可以看到这种算法比gcc的cmath库函数性能要高出20%左右。我用icc比较过，icc的速度比这个快了4倍，不知道icc具体怎么实现的。不过我推测可能是用了SIMD，因为我主要目的还是利用SIMD计算多个exp函数值，因此这里就没有考虑利用SIMD优化单个exp的计算，不知道利用SIMD能不能接近icc的性能。

log函数(ln)的实现

同样使用泰勒展开，由于$ln(0)$ln(0)的是$-\infty$−∞因此不能在0处展开，选在1处展开
$ln(1+x)=\frac{1}{n}x-\frac{1}{n^2}x^2+\frac{1}{n^3}x^3-....$ln(1+x)=n1​x−n21​x2+n31​x3−....
同样考虑两个问题，级数的收敛条件和截断误差。

$-1<x<=1$−1<x<=1级数收敛，于是我们会发现$x>1$x>1即对$ln(a),a>2$ln(a),a>2时级数是不收敛的，这时泰勒展开是得不得正确结果的。不过这并不是什么问题，因为我们还有公式$ln(a)=-ln(\frac{1}{a})$ln(a)=−ln(a1​)这个公式可以把大于2的数变成小于1的数，这时级数就收敛了。但是截断误差就比较麻烦了，以$ln(2)$ln(2)为例，由于泰勒展开中分母不再是阶乘的形式，所以这个级数收敛的异常缓慢，其误差是$O(\frac{x^2}{n})$O(nx2​)，x=1时，如果要达到$10^{-14}$10−14的精度，需要计算$10^{14}$1014项(其实是$10^7$107项)，这是不可接受的。

于是我们同样考虑exp一样的方式，把x分解为两部分$x=i*m$x=i∗m这里i满足$ln(i)$ln(i)为整数，m满足$1<=m<e$1<=m<e。这时$ln(x)=ln(i)+ln(m)$ln(x)=ln(i)+ln(m)表面上看好像并不复杂，但是实际上这并不容易实现。但是根据浮点数的特性如果i满足$log_2(i)$log2​(i)是整数，我们很容易得到$log_2(i)$log2​(i)的值，就是其指数部分，而把其指数归零即可得到对应的落在$[1,2)$[1,2)上的m(注意浮点数的指数并非正常的整型表达)。而利用换底公式，我们又可以很容易由$log_2(i)$log2​(i)得到$ln(i)$ln(i)。于是我们先进行如下操作

int _pre_log_e(double x)
{
	long unsigned int *j=((long unsigned int*)(&x));
	(*j)=(*j)<<1;
	(*j)=(*j)>>53;
	return *j-1023;
}
double _pre_log_m(double x)
{
	long unsigned int *k=((long unsigned int*)(&x));
	*k=(*k)&(0x8fffffffffffffff);
	*k=(*k)|(0x3ff0000000000000);
	return x;
}

这样我们就得到了$log_2(i)$log2​(i)和m的值，这样得到的m满足$1<=m<2$1<=m<2，$log_2(i)$log2​(i)部分已经计算完成，接下来就是计算$ln(m)$ln(m)了。虽然这时m不会等于2了，但是观察泰勒展开，可以发现当m接近2时结果同样很难快速收敛。于是我们再仿照exp的做法对m开平方，开一次平方可以把m的上限减小到1.4左右，这时$ln(m)=2ln(\sqrt{m})$ln(m)=2ln(m​)，大概最多计算40项即可达到$10^{-14}$10−14的精度，而再开方可以再减小计算的项数。由于开方运算本身有硬件支持，所以相对而言性能是比较高的，因此通过适当开方后可以显著减小计算量。根据前面的泰勒展开容易得到

$ln(\frac{1+x}{1-x})=2x+\frac{2}{3}x^3+...+\frac{2}{2n+1}x^{2n+1}$ln(1−x1+x​)=2x+32​x3+...+2n+12​x2n+1

令$m=\frac{1+x}{1-x}$m=1−x1+x​可以得到$x=\frac{m-1}{m+1}$x=m+1m−1​，这时即便m接近1.4只需要大约达到16阶即可获得$10^{-14}$10−14的精度，这时仅需计算8项。于是结合上面的结果，得到如下代码

const double _ln2=0.6931471805599453;
double _log1_2(const double &x)
{
	double dx=(x-1)/(x+1),dx2=dx*dx,r=0;
	for(int i=8;i>0;i--)
	{
		r+=1.0/(2*i+1);
		r*=dx2;
	}
	r+=1;
	r=r*2*dx;
	return r;
}
double _log(const double &x)
{
	double m=_pre_log_m(x);
	int e=_pre_log_e(x);
	m=sqrt(m);
	return 2*_log1_2(m)+e*_ln2;
}

同样与cmath的库函数比较一下

int main()
{
	double i=double(clock()%1024)/8;
	cout<<"error = "<<_log(i)-log(i)<<endl;
	double dt;
	double res[1000000];
	clock_t s=clock();	
	for(int i=1;i<1000000;i++)
	{
		res[i]=log(double(i)/50000);
	}
	dt=-double(s-clock())/CLOCKS_PER_SEC;
	cout<<res[7777]<<endl;
	s=clock();
	for(int i=1;i<1000000;i++)
	{
		res[i]=_log(double(i)/50000);
	}
	cout<<"_log "<<-double(s-clock())/CLOCKS_PER_SEC<<endl;
	cout<<"log "<<dt<<endl;
	cout<<res[7777]<<endl;
	
	return 0;
}

g++ -O3得到

_log 0.02012
log 0.028486

我们的计算同样比库函数更快。而且这里我发现一个有意思的现象，如果减小_log1_2函数中循环的次数，速度不仅不会变快，甚至还可能变慢，如果当i=8改为i=4时速度会显著降低。这个问题在关闭优化，打开调试时不会出现，但是一旦开启O3优化就会出现这个问题，可能是gcc的bug？