引言
我們知道，對數函式 ln(x) 可以展開為泰勒級數：
但是下面這個泰勒級數展開式收斂得更快：
經過簡單計算可知上式中 y = (x - 1) / (x + 1) 。
實現該演算法的 c# 程式
根據上面的第二個泰勒級數展開式，我們可以為 c# 的 decimal 資料型別實現如下的 log 擴充套件方法：
1 usingsystem;2
3 namespaceskyiv.extensions4 14
15 public static decimal log(this decimalx)16 24
25 static decimal logarithm(decimaly)26 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
this is free software with absolutely no warranty.
for details type `warranty'.
scale=30
define x(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
this is free software with absolutely no warranty.
for details type `warranty'.
scale=30
r=1.2217
a=0.9047
b=a*r
b1.10527199
define y(x)
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 usingsystem;2 usingskyiv.extensions;3
4 classtester5 )11 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 位有效數字，比我前幾天的「計算自然對數的快速演算法」一文介紹的演算法要好。
參考資料
自然對數e和圓周率pai
之前看過一部美劇叫做 疑犯追蹤 男主之一說過關於 當時聽到說的這些，很是贊同，但是覺得不能細想，因為，其他的無限不迴圈小數是不是也可以這樣。是求圓周的出來的數，為什麼偏偏是圓會有這樣的數，蘊含這樣的意義，直徑為1的圓，周長就是 看過另外一部電影中有句話說的是 everything comes a f...
自然對數的底數e的由來和計算方法
自然對數e最初的由來是銀行利息有關的問題。銀行存1塊錢，如果通貨嚴重膨脹的情況下，國家就會採取提高利率的方法，假定提高到100 那麼一年後的本息一共是2元。假定，膨脹更加嚴重了，銀行每半年一結息，那麼，一年後的本息就是 1 1 2 1 1 2 2.25元。再假定，如果最終銀行容易每天一結息，那麼一年...
求特殊自然數
總time limit 1000ms memory limit 65536kb description 乙個十進位制自然數,它的七進製與九進製表示都是三位數，且七進製與九進製的三位數碼表示順序正好相反。程式設計求此自然數,並輸出顯示。input 無。output 三行 第一行是此自然數的十進位制表示...