計算指數函數的算法
引言
我在上一篇隨筆中介紹了計算自然對數的高速算法。如今我們來看看計算指數函數的算法。我們知道。指數函數 ex 能夠展開為泰勒級數:
這個級數對全體實數 x 都收斂,而且在 x 接近零時收斂得比較快。
實現該算法的 C# 程序
依據前面所述的 ex 的泰勒級數展開式,能夠寫出下面 C# 程序來為 decimal 數據類型加入一個 Exp 擴展方法:
1 using System; 2 3 namespace Skyiv.Extensions 4 { 5 static class DecimalExtensions 6 { 7 static readonlydecimal expmax = 66.542129333754749704054283659m; 8 static readonly int[] mask = { 1, 2, 4, 8, 16, 32, 64 }; 9 static readonly decimal[] exps = 10 { 11 2.71828182845904523536028747135m, // exp(1) 12 7.38905609893065022723042746058m, // exp(2) 13 54.5981500331442390781102612029m, // exp(4)14 2980.95798704172827474359209945m, // exp(8) 15 8886110.52050787263676302374078m, // exp(16) 16 78962960182680.6951609780226351m, // exp(32) 17 6235149080811616882909238708.93m // exp(64) 18 }; 19 20 public static decimal Exp(this decimal x) 21 { 22 if (x > expmax) thrownew OverflowException("overflow"); 23 if (x < -66) return 0; 24 var n = (int)decimal.Round(x); 25 if (n > 66) n--; 26 decimal z = 1, y = Exponential(x - n); 27 for (int m = (n < 0) ?-n : n, i = 0; i < mask.Length; i++) 28 if ((m & mask[i]) != 0) z *= exps[i]; 29 return (n < 0) ? (y / z) : (y * z); 30 } 31 32 static decimal Exponential(decimal q) 33 { // q (almost) in [ -0.5, 0.5 ] 34 decimal y = 1, t = q; 35 for (var i = 1; t != 0; t *= q / ++i) y += t; 36 return y; 37 } 38 } 39 }
簡要說明例如以下:
- 第 7 行的 expmax 的值是 decimal.MaxValue 的自然對數的近似值,用於檢測 Exp 方法是否溢出(第 22 行)。
- 第 20 至 30 行的 Exp 擴展方法就是用來計算指數函數了。
- 該方法利用 ex+y = exey 這個公式。將參數 x 分為整數部分 n 和小數部分 x - n 來計算。
- 整數部分 n 又分解為 1、2、4、8、16、32、 64 諸數中某些的和,利用事先計算出來的常量來計算。
- 第 25 行是為了防止將 e66.5421 分解為 e67e-0.4579。這樣在計算 e67 時會溢出。而是分解為 e66e0.5421。
- 第 32 至 37 行的 Exponential 方法使用泰勒級數來計算 ex 。它的參數 q 越接近於零就計算得越快。
- 這個算法還是非常快的,第 35 行的 for 循環運行次數不會超過 22 次。
測試程序
以下就是調用 decimal 數據類型的 Exp 擴展方法的測試程序:
1 using System; 2 using Skyiv.Extensions; 3 4 class Tester 5 { 6 static void Main() 7 { 8 try 9 { 10 foreach (var x in new decimal[] { 11 -100, -66, -65, -1, 0, 1, 2.5m, 16, 66.5421m, 67 }) 12 Console.WriteLine("{0,-30}: exp({1})", x.Exp(), x); 13 } 14 catch (Exception ex) { Console.WriteLine(ex.Message); } 15 } 16 }
執行結果例如以下所看到的:
work$ dmcs Tester.cs DecimalExtensions.cs work$ mono Tester.exe 0 : exp(-100) 0.0000000000000000000000000000: exp(-66) 0.0000000000000000000000000001: exp(-65) 0.3678794411714423215955237702: exp(-1) 1 : exp(0) 2.7182818284590452353602874714: exp(1) 12.182493960703473438070175950: exp(2.5) 8886110.520507872636763023741 : exp(16) 79225838488862236701995526357 : exp(66.5421) overflow
能夠看出,在計算 e67 時發現了溢出。這是由於:
- decimal.MaxValue = 79,228,162,514,264,337,593,543,950,335
- e67 = 125,236,317,084,221,378,051,352,196,074.4365767534885274 ...
能夠看出,e67 已經超過 decimal 的最大值了。
事先計算的常數
在 DecimalExtensions.cs 程序的第 9 至 18 行中的 exps 靜態僅僅讀數組中存放了 e1、e2、e4、e8、e16、e32 和 e64 的值。這些值是怎樣得到的呢?這非常easy,Linux 操作系統中有一個高精度計算器 bc 。
我們能夠先編輯一個例如以下內容的文本文件 exps_in.txt:
scale=30 e(1) e(2) e(4) e(8) e(16) e(32) e(64) l(2^96-1) quit
上面的 e 代表 exp。l 代表 ln,296 - 1 就是 decimal.MaxValue。
然後運行下面命令:
work$ bc -l exps_in.txt > exps_out.txt
就能夠得出例如以下內容的輸出 exps_out.txt:
2.718281828459045235360287471352 7.389056098930650227230427460575 54.598150033144239078110261202860 2980.957987041728274743592099452888 8886110.520507872636763023740781450350 78962960182680.695160978022635108224219956195 6235149080811616882909238708.928469744831391846235799914388 66.542129333754749704054283659972
稍加整理,就能夠用在上述 C# 程序中了:
- 前 7 行就是 e 的各次冪。
- 最後一行就是 decimal.MaxValue 的自然對數。
計算指數函數的算法