計算指數函數的演算法

引言

我在上一篇隨筆中介紹了計算自然對數的快速演算法。現在我們來看看計算指數函數的演算法。我們知道，指數函數 ex 可以展開為泰勒級數：

這個級數對全體實數 x 都收斂，並且在 x 接近零時收斂得比較快。

實現該演算法的 C# 程式

根據前面所述的 ex 的泰勒級數展開式，可以寫出以下 C# 程式來為 decimal 資料類型添加一個 Exp 擴充方法：

 1 using System; 2  3 namespace Skyiv.Extensions 4 { 5   static class DecimalExtensions 6   { 7     static readonly decimal 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) throw new 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 }

簡要說明如下：

1. 第 7 行的 expmax 的值是 decimal.MaxValue 的自然對數的近似值，用於檢測 Exp 方法是否溢出(第 22 行)。
2. 第 20 至 30 行的 Exp 擴充方法就是用來計算指數函數了。
3. 該方法利用 ex+y = exey 這個公式，將參數 x 分為整數部分 n 和小數部分 x - n 來計算。
4. 整數部分 n 又分解為 1、2、4、8、16、32、 64 諸數中某些的和，利用事先計算出來的常量來計算。
5. 第 25 行是為了防止將 e66.5421 分解為 e67e-0.4579，這樣在計算 e67 時會溢出。而是分解為 e66e0.5421。
6. 第 32 至 37 行的 Exponential 方法使用泰勒級數來計算 ex 。它的參數 q 越接近於零就計算得越快。
7. 這個演算法還是很快的，第 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.cswork$ mono Tester.exe0                             : 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 的值。這些值是如何得到的呢？這很簡單，Linux 作業系統中有一個高精度計算機 bc 。我們可以先編輯一個如下內容的文字檔 exps_in.txt：

scale=30e(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.7182818284590452353602874713527.38905609893065022723042746057554.5981500331442390781102612028602980.9579870417282747435920994528888886110.52050787263676302374078145035078962960182680.6951609780226351082242199561956235149080811616882909238708.92846974483139184623579991438866.542129333754749704054283659972

稍加整理，就可以用在上述 C# 程式中了：

• 前 7 行就是 e 的各次冪。
• 最後一行就是 decimal.MaxValue 的自然對數。

參考資料

1. Wikipedia: Exponential function
2. Wikipedia: Taylor series
3. Linux man pages: bc - An arbitrary precision calculator language