两种计算自然对数的算法比较

引言

前一段时间，我写了两篇计算自然对数的算法的随笔，分别使用椭圆θ函数-算术几何平均法和泰勒级数展开式来计算。那么这两种算法的性能如何呢？在参考资料[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)