本文介紹使用泰勒級數展開式來計算自然對數，并給出相應的 C# 程式。
引言
我們知道，對數函數 ln(x) 可以展開為泰勒級數：
 |  |
但是下面這個泰勒級數展開式收斂得更快：
| | |
經過簡單計算可知上式中 y = (x - 1) / (x + 1) 。
實作該算法的 C# 程式
根據上面的第二個泰勒級數展開式，我們可以為 C# 的 decimal 資料類型實作如下的 Log 擴充方法：
1 using System;
2
3 namespace Skyiv.Extensions
4 {
5 static class DecimalExtensions
6 {
7 static readonly decimal ln10 = 2.3025850929940456840179914547m;
8 static readonly decimal lnr = 0.2002433314278771112016301167m;
9
10 public static decimal Log10(this decimal x)
11 {
12 return Log(x) / ln10;
13 }
14
15 public static decimal Log(this decimal x)
16 {
17 if (x <= 0) throw new ArgumentException("Must be positive");
18 int k = 0, l = 0;
19 for (; x > 1; k++) x /= 10;
20 for (; x <= 0.1m; k--) x *= 10; // ( 0.1, 1 ]
21 for (; x < 0.9047m; l--) x *= 1.2217m; // [ 0.9047, 1.10527199 )
22 return k * ln10 + l * lnr + Logarithm((x - 1) / (x + 1));
23 }
24
25 static decimal Logarithm(decimal y)
26 { // y in ( -0.05-, 0.05+ ), return ln((1+y)/(1-y))
27 decimal v = 1, y2 = y * y, t = y2, z = t / 3;
28 for (var i = 3; z != 0; z = (t *= y2) / (i += 2)) v += z;
29 return v * y * 2;
30 }
31 }
32 } 在這個程式中：
- 第 7 行是事先計算出來的 ln(10) 的值，用于第 12 行和第 22 行。
- 第 8 行是事先計算出來的 ln(1.2217) 的值，用于第 22 行。
- 第 15 至 23 行的 Log 擴充方法就是用來計算自然對數了。
- 通過第 19 至 20 行，将參數 x 的值變換到 ( 0.1, 1 ] 的區間中。這兩個循環隻會執行其中的一個，且循環次數不超過 28 次。
- 通過第 21 行，進一步将參數 x 的值變換到 [ 0.9047, 1.10527199 ) 的區間中。這個循環執行次數不超過 11 次。
- 第 22 行通過調用 Logarithm 方法來計算自然對數。傳入的參數是 (x - 1) / (x + 1)，其範圍大約在 ( -0.05, 0.05 ) 的區間中。
- 第 22 行的表達式是基于 ln(xy) = ln(x) + ln(y) 和 ln(xn) = n ln(x) 這兩條對數函數的運算規則。當然，後者是前者的特例。
- 第 25 至 30 行的 Logarithm 方法使用泰勒級數來計算自然對數，它的參數 y 越接近零收斂得越快。
- 注意，它的傳回值是 ln((1+y)/(1-y))，而不是 ln(y)。
- 這個算法還是很快的，第 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]。分析如下：
- 我們的目标是要将第 25 行的 Logarithm 方法的參數 y 控制在 ( -0.05, 0.05 ) 區間範圍内。
- 由前面引言中知道，x = (1+y) / (1-y)。是以計算出 x 大約在 ( 0.9047, 1.10526 ) 區間範圍内。
- 為了第 21 行的将 x 值變換到上述區間，計算出變換因子 1.10526 / 0.9047 ≈ 1.2217 。
- 這就得到第 8 行的 ln(1.2217) 的值。注意該值最後幾位是 ...6698，舍入到 ...6700，誤差相當小。(decimal 要求舍入到 28 個有效數字)
- 我原來在第 2 步采用區間 ( -0.90476, 1.105263 )，計算出來的變換因子是 1.105263 / 0.90476 ≈ 1.2216 。
- 相應的 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$
說明如下：
- 我們得到 x 在 [ 0.9047, 1.10527199 ) 區間範圍内。
- 由前面的引言可知，y = (x - 1) / (x + 1)，大約在 ( -0.05, 0.05 ) 區間範圍内。
- 這說明我們以前的計算是正确的。
測試程式
下面是調用 decimal 資料類型的 Log 和 Log10 擴充方法的測試程式：
1 using System;
2 using Skyiv.Extensions;
3
4 class Tester
5 {
6 static void Main()
7 {
8 foreach (var x in new decimal[] { 4 / decimal.MaxValue,
9 0.0000001m, 0.0001m, 0.1m, 1, 1.2217m, 2, 10, 10000,
10 100000000, decimal.MaxValue })
11 {
12 Console.WriteLine("x : " + x);
13 Console.WriteLine("ln: " + x.Log());
14 Console.WriteLine("lg: " + x.Log10());
15 Console.WriteLine();
16 }
17 }
18 } 運作結果如下所示：
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 位有效數字，比我前幾天的“計算自然對數的快速算法”一文介紹的算法要好。
參考資料
- Wikipedia: Natural logarithm
- Wikipedia: Taylor series
- Linux man pages: bc - An arbitrary precision caculator language
- 部落格園：計算自然對數的快速算法
- 部落格園：計算指數函數的算法