兩種計算自然對數的演算法比較

引言

前一段時間，我寫了兩篇計算自然對數的演算法的隨筆，分別使用橢圓θ函數-算術幾何平均法和泰勒級數展開式來計算。那麼這兩種演算法的效能如何呢？在參考資料[3]中有以下說法：

上面的 elliptic method 就是橢圓θ函數-算術幾何平均法，Taylor's method 2 就是我使用的泰勒級數展開式。可以看出，elliptic method 在計算精度大時占絕對優勢，但在計算精度小時並不佔優。而在我們的應用中，要計算的 decimal 資料類型的精度只有 28 位有效數字，即上面的 P = 28。看來在我們的應用中應該是 Taylor's method 2 佔優。

測試程式

那麼就讓我們寫一個 C# 程式來測試一下吧，下面就是 LogFactorialTest.cs：

 1 using System; 2 using System.Diagnostics; 3 using Skyiv.Extensions; 4  5 namespace Skyiv.Test 6 { 7   sealed class LogFactorialTest 8   { 9     static void Main(string[] args)10     {11       var n = (args.Length > 0) ? int.Parse(args[0]) : 10;12       Run(2);13       Run(n);14     }15     16     static void Run(int n)17     {18       var sw = Stopwatch.StartNew();19       var v = LogFactorial(n);20       sw.Stop();21       Console.WriteLine("{0} {1,-30}: ln({2:N0}!)", sw.Elapsed, v, n);22     }23     24     static decimal LogFactorial(int n)25     {26       var v = 0m;27       for (var i = 1; i <= n; i++) v += ((decimal)i).Log();28       return v;29     }30   }31 }

這個程式通過計算 ln(n!) 來測試兩種計算自然對數的演算法的效能。簡要說明如下：

1. 第 11 行取得使用者在命令列參數指定的 n 值（如未指定，則使用預設值 10），用以計算 ln(n!)。
2. 第 12 行先計算一次 ln(2!) 作熱身，以便 CLR 的 JIT 先將要啟動並執行演算法代碼編譯為機器碼。另外這兩種計算自然對數的演算法計算出來的 ln(2) 的值略有不同，也可以作為一個區分的標誌。
3. 第 24 至 29 行的 LogFactorial 方法就是用來計算 ln(n!) 的。它利用 ln(n!) = ln(1) + ln(2) + ... + ln(n) 來計算。
4. 本程式使用 Stopwatch 類計時。

Linux 作業系統下編譯和運行

在 Arch Linux 64-bit 作業系統的 mono 2.10.8 環境下編譯和運行，結果如下所示。其中 DecimalExtensions1.cs 是參考資料[1]中的程式，DecimalExtensions2.cs 是參考資料[2]中的程式，LogFactorialTest.cs 就是上一小節的程式。

work$ dmcs --versionMono C# compiler version 2.10.8.0work$ dmcs LogFactorialTest.cs DecimalExtensions1.cs -out:LogFactorialTest1.exework$ dmcs LogFactorialTest.cs DecimalExtensions2.cs -out:LogFactorialTest2.exework$ mono LogFactorialTest1.exe 1000000000:00:00.0022244 0.6931471805599453094172321215: ln(2!)00:01:44.7091158 151180965.48756956489842537225: ln(10,000,000!)work$ mono LogFactorialTest2.exe 1000000000:00:00.0181903 0.6931471805599453094172321225: ln(2!)00:03:54.9390478 151180965.48756956489842537227: ln(10,000,000!)work$ mono LogFactorialTest1.exe 10000000000:00:00.0022159 0.6931471805599453094172321215: ln(2!)00:17:57.6529645 1742068084.5245154532285821925: ln(100,000,000!)work$ mono LogFactorialTest2.exe 10000000000:00:00.0133964 0.6931471805599453094172321225: ln(2!)00:39:23.8652797 1742068084.5245154532285821925: ln(100,000,000!)work$ mono LogFactorialTest1.exe 100000000000:00:00.0011895 0.6931471805599453094172321215: ln(2!)03:04:05.5218954 19723265848.226982607923141006: ln(1,000,000,000!)work$ mono LogFactorialTest2.exe 100000000000:00:00.0018197 0.6931471805599453094172321225: ln(2!)05:27:32.3909935 19723265848.226982607923141007: ln(1,000,000,000!)

我們的測試程式分別使用這兩種演算法計算了 ln(107!)、ln(108!) 和 ln(109!) 的值，計算時間最長的將近五個半小時。

Windows 7 作業系統下編譯和運行

在 Windows 7 SP1 32-bit 作業系統的 .NET Framework 4.5 環境下編譯和運行：

D:\work> csc -out:LogFactorialTest1.exe LogFactorialTest.cs DecimalExtensions1.csMicrosoft(R) Visual C# 編譯器版本 4.0.30319.17929用於 Microsoft(R) .NET Framework 4.5著作權 (C) Microsoft Corporation。著作權所有，並保留一切權利。D:\work> csc -out:LogFactorialTest2.exe LogFactorialTest.cs DecimalExtensions2.csMicrosoft(R) Visual C# 編譯器版本 4.0.30319.17929用於 Microsoft(R) .NET Framework 4.5著作權 (C) Microsoft Corporation。著作權所有，並保留一切權利。D:\work> LogFactorialTest1 1000000000:00:00.0034542 0.6931471805599453094172321215: ln(2!)00:01:00.1048788 151180965.48756956489842537224: ln(10,000,000!)D:\work> LogFactorialTest2 1000000000:00:00.0043189 0.6931471805599453094172321214: ln(2!)00:01:47.1634292 151180965.48756956489842537225: ln(10,000,000!)D:\work> LogFactorialTest1 10000000000:00:00.0034569 0.6931471805599453094172321215: ln(2!)00:09:21.1743684 1742068084.5245154532285821925: ln(100,000,000!)D:\work> LogFactorialTest2 10000000000:00:00.0045334 0.6931471805599453094172321214: ln(2!)00:18:13.4201181 1742068084.5245154532285821924: ln(100,000,000!)D:\work> LogFactorialTest1 100000000000:00:00.0035446 0.6931471805599453094172321215: ln(2!)01:36:06.8523762 19723265848.226982607923141006: ln(1,000,000,000!)D:\work> LogFactorialTest2 100000000000:00:00.0043396 0.6931471805599453094172321214: ln(2!)03:05:45.9748574 19723265848.226982607923141006: ln(1,000,000,000!)

可以看出，同樣的程式，在這台機器上運行速度更快。這兩台機器的型號是一樣的，但這台機器 CPU 的主頻比前一台稍高。

運行結果分析

上面兩個運行結果整理如下表所示：

可見，在我們的應用中，演算法2(橢圓θ函數-算術幾何平均法)比演算法1(泰勒級數展開式)大約要慢一倍左右。

運行環境

第一台機器是 2010-10-13 出廠的 Lenovo ThinkCentre M6100t PC 機，軟體和硬體的資訊如下所示：

第二台機器的型號和第一台相同，但出廠期日稍遲，是 2011-12-02 出廠的。相應的資訊如下所示：

參考資料

1. 部落格園：計算自然對數的演算法
2. 部落格園：計算自然對數的快速演算法
3. CiNii Articles: Practically Fast Multiple-Precision Evaluation of LOG(X)