使用公式
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; } }
留言列表