实现log()和exp()函数的方法,并以此计算pow()

由于本人水平极度有限，所以所写的只是个人的想法，不可避免地存在错漏，如果有任何意见或批评，非常欢迎告诉我，请email到pjy.studio@gmail.com

演示代码下载http://pjy.studio.googlepages.com/log_exp_pow.cpp      

如果不能正常查看文章，请点击下载文章

最近特然对如何自己写出log()函数十分感兴趣，于是就想探索一下其实现方法。

上次简单的分析了一下pow()函数的实现，在那里只是简单的运用了Taylor展开。鉴于自然对数函数log(),函数exp()和函数pow()有密切的关系,于是就有了一齐把这三个函数的实现方法都写出来的想法。而其中，实现log()函数是最关健和复杂的部分。为了与库函数分开，我自己写的这三个函数名命名为powP，lnP，expP。以下是我的一些想法：

    一，运用Taylor展开

要实现lnP，应用Taylor展开是最一般的方法了。这里把lnP(t)转化为lnP(1+x)，然后进行展开，由数学分析知识知道在0<t<=2的范围类，其展开式是收敛的，记 。这样就可以粗略计算lnP了，但实现时会发现，用Taylor展开计算的效果并不令人满意，这主要是指其收敛的速度不够快，当t接近0，1，2时，计算时间十分长。这样，我们就需要对   做出变换，以加快其收敛速度，这里，我用了Aitken外推数列的方法。

 

二，应用Aitken算法加速数列收敛

在实现lnP函数中，加速数列收敛的Aitken方法是整个过程的核心，因此，简单的介绍一下Aitken方法：

计算级数 的和 a ，记级数部分和为  ，如果此数列收敛慢，则有必要构造新数列

，使他更快地收敛到 a 。记  ，若   的收敛阶为1，即   ， c 为常数，于是

 

，解之有    ，于是，便得到 Aitken 外推数列

 

 

，记   ，要证明外推数列比原数列更快的收敛于 a ，只需满足

 

 

及   即可。可以证明 lnP 的 Taylor 展开满足要求。（由于书写不便，证明不给出了，关于 Aitken 方法，参考的是《数值试验》一书）。这样只要应用 Aitken 外推出新的数列，计算 lnP 的速度就大大加快。

 

 三，进一步改进

虽然应用Aitken方法，得出收敛更快的数列，但通过试验，可以发现当t十分接近0时，其收敛速度依然不能令人满意。而这是个十分重要的问题，因为应用基于Taylor展开的方法，只能计算0<t<=2的lnP(t).而当t大于2时，就要转换为计算-lnP(1/t)，这就把问题转为计算0<t<0.5的lnP(t)，如果此时不能解决t接近0时收敛慢的问题，也就相当于解决不了t比较大时的计算问题。

  要解决这个问题，可以只需把计算在(0,0.5)的t转为计算在(1,2)的t。要实现这样的转换，可以进行t=a×e^n的转化，其中a在(1,2)内。这样，lnP(t)=ln(a)+n，只要找出对应的n即可。对于小于0.5的t，要找到满足1<t/e^n<2的n，只要定义一个n的初值，然后对n不断减少，计算n为何值时，t/e^n=a在(1,2)即可。为了确保每次减少n，得出的两个范围（e^n,2*e^n）始终不会分离，我们要解2*e^(n-d) > e^n。可以计算得，d取得比较大的值为0.6。所以可以取n= -1，这时满足1<t/e^n<2的t范围是（0.3678794，0.7357588），例如每次减少0.6，当减少一次时，n= -1.6，此时满足1<t/e^n<2的t范围是（0.2018965，0.403793），减少第二次时，n= -2.2，范围为（0.1108031，0.2216063）。这样，只要如此计算，就能找到合适的n，使1<t/e^n<2。这样的话，又需要有计算e^n的函数，我们可以用相同的方法，先进行Taylor展开，在用Aitken外推数列，即可得出计算e^n的函数expP()。

综上述，已能较好的解决t在(0,2)内的自然对数，最后对于t>2进行-lnP(1/t)转换，就能计算t>0的自然对数。

最后，有了lnP()与expP()，就可以计算powP()了，因为ln(x^y) = y×lnx ,所以x^y=exp(y×lnx )。

由以上分析，可以写出一下代码：

 

 

double  expP( double  x) // 计算e^x,实现系统的exp()功能 
... {
    if(x==0) return 1;
    if(x<0) return 1/expP(-x); 
    double y=x,ex_p1=0,ex_p2=0,ex_p3=0,ex_p=0,ex_px=0,ex_tmp=1,dex_px=1,tmp;
    int l;
    for(l=1,tmp=1;((ex_px-ex_tmp)>1e-10 || (ex_px-ex_tmp)<-1e-10) && dex_px>1e-10;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  lnP( double  x)
... {
    if(x==1) return 0;
    else if(x>2) return -lnP(1/x);
    else if(x<.1)
    ...{
     double n=-1;
     double a;
     do
     ...{
     n=n-.6;
     a=x/exp(n);
     }
     while(a>2 || a<1);
     return lnP(a)+n;
    }
    double y=x-1,ln_p1=0,ln_p2=0,ln_p3=0,ln_p=0,ln_px=0,ln_tmp=1,dln_px=1,tmp;
    int l;
    for(l=1,tmp=1;(ln_px-ln_tmp)>1e-10 || (ln_px-ln_tmp)<-1e-10;l++)
    ...{
        ln_tmp=ln_px;
        tmp*=y;
        if(l==1) tmp=tmp/l;
        else tmp=tmp/-l;
        ln_p1+=tmp;
        ln_p2=ln_p1+-1*tmp*y*l/(l+1);
        ln_p3=ln_p2+tmp*y*y*l/(l+2);
        dln_px=ln_p3-ln_p2;
        ln_px=ln_p3-dln_px*dln_px/(ln_p3-2*ln_p2+ln_p1);
        tmp*=l;
        }
        return ln_px;
}
double  powP( double  x, double  y) // 计算x^y 
... {
     if(x==0 && y!=0) return 0;
     else if(x==0 && y==0) return 1;
     else if(y<0) return 1/powP(x,-y);//把指数小于0的情况转为1/x^-y计算
     else if(x<0 && y-int(y)!=0) return 0;//若x为负，且y不为整数数，则出错，返回0  
     else if(x<0 && y-int(y)==0)//若x为负，且y为整数数，则用循环计算 
          ...{
          double powint=1;
          int i;
          for(i=1;i<=y;i++) powint*=x;
          return powint;
          }
     return expP(y*lnP(x));
}