Edison.X. Blog

使用公式

exp(x) = 1 + x + x^2 / 2! + x^3 / 3! + ... + x^n / n!

不要用定義式 : exp(x) = (1+x/n)^n , n->0

誤差沒辦法控制。

誤差作收斂條件

純粹照著公式跑。

double xExp(const double x, const double eps)
{
     double ans1 ;
     double ans2 = 1.0;
     double fact=1, xn=x, cnt=2.0;
     do{
           ans1=ans2, ans2 = ans1+xn/fact;
           fact*=cnt;
           xn = xn*x;
           cnt+=1.0;
     }while( (ans1 > ans2 + eps) || (ans2 > ans1 + eps));
     return ans2;
}

固定次數為收斂條件

將公式改成

exp(x) = 1 + x + x^2 / 2! + x^3 / 3! + ... + x^n / n!

= 1 + x * (1 + x/2 * (1+ x/3 * (1+x/4 * ....*(1+x/n)...)))

由於必須從內層算到外層，故須先確定 n 有幾項，沒辦法使用誤差作收斂條件。

double xExpn(const double x, const int itera)
{
     double ret;    
     int n=itera;

     ret=1.0+x/n, --n;
     while(n) ret = 1.0+ret*x/n, --n;
     return ret;
}

Fast Exp

在 x 數值很小、允許約有 5% 左右誤差時，將 taylor series 只取前幾項出來，並將分母通分提出。

inline float fastexp3(float x) {return (6+x*(6+x*(3+x)))*0.16666666f;}

inline float fastexp4(float x) {return (24+x*(24+x*(12+x*(4+x))))*0.041666666f;}

inline float fastexp5(float x) {return (120+x*(120+x*(60+x*(20+x*(5+x)))))*0.0083333333f;}

inline float fastexp6(float x) {return (720+x*(720+x*(360+x*(120+x*(30+x*(6+x))))))*0.0013888888f;}

inline float fastexp7(float x) {return (5040+x*(5040+x*(2520+x*(840+x*(210+x*(42+x*(7+x)))))))*0.00019841269f;}

inline float fastexp8(float x) {return (40320+x*(40320+x*(20160+x*(6720+x*(1680+x*(336+x*(56+x*(8+x))))))))*2.4801587301e-5;}

被吃掉的話看這份

inline float fastexp3(float x) {return (6+x*(6+x*(3+x)))*0.16666666f;}

inline float fastexp4(float x) {return (24+x*(24+x*(12+x*(4+x))))*0.041666666f;}

inline float fastexp5(float x) {return (120+x*(120+x*(60+x*(20+x*(5+x)))))*0.0083333333f;}

inline float fastexp6(float x) {return (720+x*(720+x*(360+x*(120+x*(30+x*(6+x))))))*0.0013888888f;}

inline float fastexp7(float x) {return (5040+x*(5040+x*(2520+x*(840+x*(210+x*(42+x*(7+x)))))))*0.00019841269f;}

inline float fastexp8(float x) {return (40320+x*(40320+x*(20160+x*(6720+x*(1680+x*(336+x*(56+x*(8+x))))))))*2.4801587301e-5;}

此處採用 fastexp9 部份

double xExpSamll(const double x)
{
     return (362880+x*(362880+x*(181440+x*(60480+x*(15120+x*(3024+x*(504+x*(72+x*(9+x)))))))))*2.755731922398589e-006;
}

截至目前，上面的所有程式碼都不算實用。

考慮實際因素

1. 當 x 愈大時，收斂速度慢。

2. 實際上和 math.h 裡面之 exp ，能表達數值少很多。

鑑於 x 在值小時，收斂速度較快 (計算次數較少)，試著將 x 值化開

令 x = I + F，其中 I 為整數部份，F 為浮點數部份，exp(I+F) = exp(I) * exp(F)

而 exp(I) 部份可調用 fastpower ，只需 O( log(I) ) 即可完成，

F 部份則再進入上述之 xExp 即可，最後將結果相乘，可表示範圍也與 exp 相差無幾。

有個動作建議事先判斷，當 x > 709 時，exp 會暴，

x < -709 時會收到 0 ，便不會在 loop 裡面空耗時間。

// ------------------------------------------------------------
double xFastPow(const double base, const int exponment)
{
     double rst=1.0, pwr=base;
     int e = (exponment > 0) ? (exponment) : (-exponment);
     while(e){
           if(e&1) rst*=pwr;
           pwr*=pwr;
           e>>=1;
     }
     return (exponment > 0) ? rst : 1.0/rst;
}
// ------------------------------------------------------------
double xExpFast(const double exponment, double eps)
{
     const double Euler = 2.718281828459045;
     double rst=1.0;
     double p1=ceil(exponment), p2=exponment-p1;

     if(exponment > 709.0){
           p1=1.0, p2=0.0;
           return p1/p2;
     } else if( exponment < -709.0) {
           return 0.0;
     } else {
           return xExp(p2, eps) * xFastPow(Euler, p1) ;
     }   
}

此法效率約 compiler (gcc / vc) 慢 2~5 倍，但上述效率慢 50 倍以上。

再改善 ?

硬要再改善的方式有兩種

(1) 既然已化成 exp(I + F)，若對於精度不是非常要求，可在 exp(F) 部份調用 xExpSmall 加快速度。

(2) 既然 exp(I) 都已確定只有 709 個數字，可考慮建 table 0~709 ，這部份變成 O(1)。

參考如下

double xExpTable(const double exponment)
{
     const double Euler = 2.718281828459045;
     double rst=1.0;
     double p1=ceil(exponment), p2=exponment-p1;

     if(exponment > 709.0){
           p1=1.0, p2=0.0;
           return p1/p2;
     } else if( exponment < -709.0) {
           return 0.0;
     } else {
           rst = (p1 < 0) ? (1.0 / exp_table[-(int)p1]-1) : (exp_table[(int)p1]);
           return xExpSamll(p2) * rst;
     }   
}