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

引言

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

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

測試程序

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

using System;
using System.Diagnostics;
using Skyiv.Extensions;

namespace Skyiv.Test
{
  sealed class LogFactorialTest
  {
    static void Main(string[] args)
    {
      var n = (args.Length > 0) ? int.Parse(args[0]) : 10;
      Run(2);
      Run(n);
    }
    
    static void Run(int n)
    {
      var sw = Stopwatch.StartNew();
      var v = LogFactorial(n);
      sw.Stop();
      Console.WriteLine("{0} {1,-30}: ln({2:N0}!)", sw.Elapsed, v, n);
    }
    
    static decimal LogFactorial(int n)
    {
      var v = 0m;
      for (var i = 1; i <= n; i++) v += ((decimal)i).Log();
      return v;
    }
  }
}

這個程序通過計算 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 --version
Mono C# compiler version 2.10.8.0
work$ dmcs LogFactorialTest.cs DecimalExtensions1.cs -out:LogFactorialTest1.exe
work$ dmcs LogFactorialTest.cs DecimalExtensions2.cs -out:LogFactorialTest2.exe
work$ mono LogFactorialTest1.exe 10000000
00:00:00.0022244 0.6931471805599453094172321215: ln(2!)
00:01:44.7091158 151180965.48756956489842537225: ln(10,000,000!)
work$ mono LogFactorialTest2.exe 10000000
00:00:00.0181903 0.6931471805599453094172321225: ln(2!)
00:03:54.9390478 151180965.48756956489842537227: ln(10,000,000!)
work$ mono LogFactorialTest1.exe 100000000
00:00:00.0022159 0.6931471805599453094172321215: ln(2!)
00:17:57.6529645 1742068084.5245154532285821925: ln(100,000,000!)
work$ mono LogFactorialTest2.exe 100000000
00:00:00.0133964 0.6931471805599453094172321225: ln(2!)
00:39:23.8652797 1742068084.5245154532285821925: ln(100,000,000!)
work$ mono LogFactorialTest1.exe 1000000000
00:00:00.0011895 0.6931471805599453094172321215: ln(2!)
03:04:05.5218954 19723265848.226982607923141006: ln(1,000,000,000!)
work$ mono LogFactorialTest2.exe 1000000000
00:00:00.0018197 0.6931471805599453094172321225: ln(2!)
05:27:32.3909935 19723265848.226982607923141007: ln(1,000,000,000!)

我們的測試程序分別使用這兩種算法計算了 ln(10⁷!)、ln(10⁸!) 和 ln(10⁹!) 的值，計算時間最長的將近五個半小時。

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

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

D:\work> csc -out:LogFactorialTest1.exe LogFactorialTest.cs DecimalExtensions1.cs
Microsoft(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.cs
Microsoft(R) Visual C# 編譯器版本 4.0.30319.17929
用於 Microsoft(R) .NET Framework 4.5
版權所有 (C) Microsoft Corporation。保留所有權利。
D:\work> LogFactorialTest1 10000000
00:00:00.0034542 0.6931471805599453094172321215: ln(2!)
00:01:00.1048788 151180965.48756956489842537224: ln(10,000,000!)
D:\work> LogFactorialTest2 10000000
00:00:00.0043189 0.6931471805599453094172321214: ln(2!)
00:01:47.1634292 151180965.48756956489842537225: ln(10,000,000!)
D:\work> LogFactorialTest1 100000000
00:00:00.0034569 0.6931471805599453094172321215: ln(2!)
00:09:21.1743684 1742068084.5245154532285821925: ln(100,000,000!)
D:\work> LogFactorialTest2 100000000
00:00:00.0045334 0.6931471805599453094172321214: ln(2!)
00:18:13.4201181 1742068084.5245154532285821924: ln(100,000,000!)
D:\work> LogFactorialTest1 1000000000
00:00:00.0035446 0.6931471805599453094172321215: ln(2!)
01:36:06.8523762 19723265848.226982607923141006: ln(1,000,000,000!)
D:\work> LogFactorialTest2 1000000000
00: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)