引言
在1982年,Tateaki. Sasaki 和 Yasumasa Kanada 發表了一篇論文:Practically Fast Multiple-Precision Evaluation of LOG(x)。在這篇只有四頁的論文中,他們介紹了一個計算自然對數的快速算法。
實現該算法的 C# 程序
我們知道,.NET Framework Class Library 中的 System.Math.Log 方法只能對 dobule 數據類型計算自然對數。現在我們為 decimal 數據類型實現一個 Log 擴展方法,如下所示:
1 using System; 2 3 namespace Skyiv.Extensions 4 { 5 static class DecimalExtensions 6 { 7 static readonly decimal pi2 = 6.283185307179586476925286766559m; 8 static readonly decimal ln10 = 2.30258509299404568401799145468m; 9 static readonly decimal eps2 = 0.00000000000001m; // 1e-14 10 static readonly decimal eps1 = eps2 * eps2; // 1e-28 11 12 public static decimal Sqrt(this decimal x) 13 { 14 if (x < 0) throw new ArgumentException("Must be nonnegative"); 15 if (x == 0) return 0; 16 var y = (decimal)Math.Sqrt((double)x); 17 return (y + x / y) / 2; 18 } 19 20 public static decimal Log10(this decimal x) 21 { 22 return Log(x) / ln10; 23 } 24 25 public static decimal Log(this decimal x) 26 { 27 if (x <= 0) throw new ArgumentException("Must be positive"); 28 if (x == 1) return 0; 29 var k = 0; 30 for (; x > 0.1m; k++) x /= 10; 31 for (; x <= 0.01m; k--) x *= 10; 32 return k * ln10 - NegativeLog(x); 33 } 34 35 static decimal NegativeLog(decimal q) 36 { // q in ( 0.01, 0.1 ] 37 decimal r = q, s = q, n = q, q2 = q * q, q1 = q2 * q; 38 for (var p = true; (n *= q1) > eps1; s += n, q1 *= q2) 39 r += (p = !p) ? n : -n; 40 decimal u = 1 - 2 * r, v = 1 + 2 * s, t = u / v; 41 decimal a = 1, b = Sqrt(1 - t * t * t * t); 42 for (; a - b > eps2; b = Sqrt(a * b), a = t) t = (a + b) / 2; 43 return pi2 / (a + b) / v / v; 44 } 45 } 46 }
在上述程序中:
- 第 12 至 18 行實現 Sqrt 擴展方法。使用牛頓迭代法,僅迭代一次 (第 17 行)。這是因為第 16 行使用 Math.Sqrt (double) 方法給出的初值的精度已經達到 15,而 decimal 的精度為 28,迭代一足夠了。
- 第 25 至 33 行實現 Log 擴展方法。該方法先將參數 x 變換到 ( 0.01, 0.1 ] 區間中 ( 第 29 至 31 行),然后調用 NegativeLog 方法進行計算。
- 第 35 至 44 行的 NegativeLog 方法就是我們算法的核心,使用 橢圓θ函數 ( 第 37 至 40 行 ) 和 算術幾何平均法 ( 第 41 至 43 行 ) 來快速計算自然對數。該算法的原理請參閱參考資料[1]。
- 第 7 行是事先計算出來的 2π 的值,在第 43 行中使用。
- 第 8 行是事先計算出來的 ln10 的值,在第 32 行和 22 行中使用。
- 我們的程序中使用 ln10 的值,而參考資料[1]中使用 ln2 的值。這是因為 decimal 使用十進制比較好。而一般的 double 使用二進制比較好。
- 第 10 行和第 9 行的 eps1 = 10-28 和 eps2 = 10-14 使用十進制控制計算精度。而參考資料[1]中使用二進制控制計算精度,即 2-p 和 2-p/2。
算法的性能
在上述程序中加入一些調試語句,就可以在程序運行時報告一些重要變量的值,結果如下所示:
q : 0.1000000000000000000000000000 k: 2 r : 0.0999000009999999000000001000 s : 0.1001000010000001000000001000 u : 0.8001999980000001999999998000 v : 1.2002000020000002000000002000 loop1: 4 b0: 0.8957696680606997488130397041 a : 0.9471678269798758062616205879 b : 0.9471678269798660936897543331 loop2: 3 NL: 2.3025850929940456840179914544 Ln: 2.3025850929940456840179914550 q : 0.0100000000000000000000000001 k: 28 r : 0.0099999900000000010000000001 s : 0.0100000100000000010000000001 u : 0.9800000199999999979999999998 v : 1.0200000200000000020000000002 loop1: 2 b0: 0.3845443947629648387969403232 a : 0.6556979528724023734005530455 b : 0.6556979528723956496570849678 loop2: 4 NL: 4.6051701859880913680359828988 Ln: 59.867212417845187784467777833
上面是對 10 和 100000000000000000000000001 分別調用 Log 方法計算自然對數的調試結果:
- k 是 Log 方法執行到第 32 行時的值。其他變量都是第 35 至 44 行中 NegativeLog 方法執行時的值。
- loop1 是第 38 至 39 行的 for 循環的執行次數。
- loop2 是第 42 行的 for 循環執行次數。
- q, r, s, u, v, a, b 都是計算終了時的值。
- b0 是 b 的初值 (第 41 行)。a 的初值是 1 。
- NL 是 NegativeLog 方法的返回值,即參數 q 的自然對數的相反數。
- Ln 是 Log 方法的返回值,即參數 x 的自然對數。
- 從上述調試結果可以看出,這個算法是非常高效的。算法中的兩個 for 循環執行次數一般都不超過 4。
測試程序
下面是調用 decimal 數據類型的 Sqrt、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.1m, 1, 10, 100000000, decimal.MaxValue }) 10 { 11 Console.WriteLine("x : " + x); 12 Console.WriteLine("sqrt: " + x.Sqrt()); 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 sqrt: 0.00000000000001 ln : -64.472382603833279152503760731 lg : -28.000000000000000000000000000 x : 0.0000001 sqrt: 0.0003162277660168379331998894 ln : -16.118095650958319788125940182 lg : -6.9999999999999999999999999996 x : 0.1 sqrt: 0.3162277660168379331998893545 ln : -2.3025850929940456840179914544 lg : -0.9999999999999999999999999999 x : 1 sqrt: 1 ln : 0 lg : 0 x : 10 sqrt: 3.1622776601683793319988935445 ln : 2.3025850929940456840179914550 lg : 1.0000000000000000000000000001 x : 100000000 sqrt: 10000 ln : 18.420680743952365472143931638 lg : 8.000000000000000000000000000 x : 79228162514264337593543950335 sqrt: 281474976710656.00000000000000 ln : 66.542129333754749704054283660 lg : 28.898879583742194740518933893
從中可以看出,這個算法的運算精度能夠達到 27,只有最后一位可能有誤差。