計算自然對數的算法

引言

我們知道，對數函數 ln(x) 可以展開為泰勒級數：

但是下面這個泰勒級數展開式收斂得更快：

經過簡單計算可知上式中 y = (x - 1) / (x + 1) 。

實現該算法的 C# 程序

根據上面的第二個泰勒級數展開式，我們可以為 C# 的 decimal 數據類型實現如下的 Log 擴展方法：

using System;

namespace Skyiv.Extensions
{
  static class DecimalExtensions
  {
    static readonly decimal ln10 = 2.3025850929940456840179914547m;
    static readonly decimal lnr = 0.2002433314278771112016301167m;

    public static decimal Log10(this decimal x)
    {
      return Log(x) / ln10;
    }

    public static decimal Log(this decimal x)
    {
      if (x <= 0) throw new ArgumentException("Must be positive");
      int k = 0, l = 0;
      for (; x > 1; k++) x /= 10;
      for (; x <= 0.1m; k--) x *= 10;        // ( 0.1, 1 ]
      for (; x < 0.9047m; l--) x *= 1.2217m; // [ 0.9047, 1.10527199 )
      return k * ln10 + l * lnr + Logarithm((x - 1) / (x + 1));
    }
    
    static decimal Logarithm(decimal y)
    { // y in ( -0.05-, 0.05+ ), return ln((1+y)/(1-y))
      decimal v = 1, y2 = y * y, t = y2, z = t / 3;
      for (var i = 3; z != 0; z = (t *= y2) / (i += 2)) v += z;
      return v * y * 2;
    }
  }
}

在這個程序中：

1. 第 7 行是事先計算出來的 ln(10) 的值，用於第 12 行和第 22 行。
2. 第 8 行是事先計算出來的 ln(1.2217) 的值，用於第 22 行。
3. 第 15 至 23 行的 Log 擴展方法就是用來計算自然對數了。
4. 通過第 19 至 20 行，將參數 x 的值變換到 ( 0.1, 1 ] 的區間中。這兩個循環只會執行其中的一個，且循環次數不超過 28 次。
5. 通過第 21 行，進一步將參數 x 的值變換到 [ 0.9047, 1.10527199 ) 的區間中。這個循環執行次數不超過 11 次。
6. 第 22 行通過調用 Logarithm 方法來計算自然對數。傳入的參數是 (x - 1) / (x + 1)，其范圍大約在 ( -0.05, 0.05 ) 的區間中。
7. 第 22 行的表達式是基於 ln(xy) = ln(x) + ln(y) 和 ln(xⁿ) = n ln(x) 這兩條對數函數的運算規則。當然，后者是前者的特例。
8. 第 25 至 30 行的 Logarithm 方法使用泰勒級數來計算自然對數，它的參數 y 越接近零收斂得越快。
9. 注意，它的返回值是 ln((1+y)/(1-y))，而不是 ln(y)。
10. 這個算法還是很快的，第 28 行的 for 循環執行次數不會超過 10 次。

程序中相關常數的由來

上面程序中的 1.2217 和 0.9047 等常數是如何得到的呢？請看下面的計算：

work$ bc -l
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
scale=30
define x(y) { return (1+y)/(1-y); }
x(-0.05)
.904761904761904761904761904761
x(0.05)
1.105263157894736842105263157894
1.10526/0.9047
1.221686746987951807228915662650
l(1.2217)
.200243331427877111201630116698
l(1.2216)
.200161474922285626409839638619
quit
work$

上面使用 Linux 中的 bc 進行計算，l 代表 ln 函數，請參閱參數資料[3]。分析如下：

1. 我們的目標是要將第 25 行的 Logarithm 方法的參數 y 控制在 ( -0.05, 0.05 ) 區間范圍內。
2. 由前面引言中知道，x = (1+y) / (1-y)。所以計算出 x 大約在 ( 0.9047, 1.10526 ) 區間范圍內。
3. 為了第 21 行的將 x 值變換到上述區間，計算出變換因子 1.10526 / 0.9047 ≈ 1.2217 。
4. 這就得到第 8 行的 ln(1.2217) 的值。注意該值最后幾位是 ...6698，舍入到 ...6700，誤差相當小。(decimal 要求舍入到 28 個有效數字)
5. 我原來在第 2 步采用區間 ( -0.90476, 1.105263 )，計算出來的變換因子是 1.105263 / 0.90476 ≈ 1.2216 。
6. 相應的 ln(1.2216) 的最后幾位是 ...8619，舍入到 ...8600，誤差就稍微大了一點。

驗證常數的值

讓我們來驗證一下前面計算的常數的值是否正確：

work$ bc -l
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
scale=30
r=1.2217
a=0.9047
b=a*r
b
1.10527199
define y(x) { return (x-1)/(x+1); }
y(a)
-.050034126109098545702735338898
y(b)
.050003985470779953710399196447
quit
work$ 

說明如下：

1. 我們得到 x 在 [ 0.9047, 1.10527199 ) 區間范圍內。
2. 由前面的引言可知，y = (x - 1) / (x + 1)，大約在 ( -0.05, 0.05 ) 區間范圍內。
3. 這說明我們以前的計算是正確的。

測試程序

下面是調用 decimal 數據類型的 Log 和 Log10 擴展方法的測試程序：

using System;
using Skyiv.Extensions;

class Tester
{
  static void Main()
  {
    foreach (var x in new decimal[] { 4 / decimal.MaxValue,
      0.0000001m, 0.0001m, 0.1m, 1, 1.2217m, 2, 10, 10000,
      100000000, decimal.MaxValue })
    {
      Console.WriteLine("x : " + x);
      Console.WriteLine("ln: " + x.Log());
      Console.WriteLine("lg: " + x.Log10());
      Console.WriteLine();
    }
  }
}

運行結果如下所示：

work$ dmcs Tester.cs DecimalExtensions.cs
work$ mono Tester.exe
x : 0.0000000000000000000000000001
ln: -64.472382603833279152503760732
lg: -28.000000000000000000000000000

x : 0.0000001
ln: -16.118095650958319788125940183
lg: -7.0000000000000000000000000000

x : 0.0001
ln: -9.210340371976182736071965819
lg: -4.0000000000000000000000000001

x : 0.1
ln: -2.3025850929940456840179914547
lg: -1

x : 1
ln: 0
lg: 0

x : 1.2217
ln: 0.2002433314278771112016301167
lg: 0.0869645738770510340282719812

x : 2
ln: 0.6931471805599453094172321215
lg: 0.3010299956639811952137388947

x : 10
ln: 2.3025850929940456840179914547
lg: 1

x : 10000
ln: 9.210340371976182736071965819
lg: 4.0000000000000000000000000001

x : 100000000
ln: 18.420680743952365472143931638
lg: 8.000000000000000000000000000

x : 79228162514264337593543950335
ln: 66.542129333754749704054283660
lg: 28.898879583742194740518933893

從上面運行結果可以看出，精度基本上達到了 28 位有效數字，比我前幾天的“計算自然對數的快速算法”一文介紹的算法要好。

參考資料

1. Wikipedia: Natural logarithm
2. Wikipedia: Taylor series
3. Linux man pages: bc - An arbitrary precision caculator language
4. 博客園：計算自然對數的快速算法
5. 博客園：計算指數函數的算法