計算冪函數的幾種方法

引言

我們知道，自然對數的底 e 定義為以下極限值：

這個公式很適合於對冪函數的計算進行一些測試，得到的結果是 e 的近似值，不用擔心當 n 很大時計算結果會溢出。

測試程序

下面就是 Tester.cs：

using System;
using System.Numerics;
using System.Diagnostics;
using Skyiv.Extensions;

namespace Skyiv.Test
{
  sealed class Tester
  {
    string Standard(long n)
    { // n == 10^m
      if (n > 100000) return "Skip";
      var s = BigInteger.Pow(n + 1, (int)n).ToString();
      s = s.Substring(0, Math.Min(31, s.Length));
      return s[0] + "." + s.Substring(1);
    }
    
    string Direct(long n)
    {
      if (n > 1000000000) return "Skip";
      var y = 1m;
      for (var x = 1 + 1m / n; n > 0; n--) y *= x;
      return y.ToString();
    }
    
    string Binary(long n)
    {
      var y = 1m;
      for (var x = 1 + 1m / n; n != 0; x *= x, n >>= 1)
        if ((n & 1) != 0) y *= x;
      return y.ToString();
    }
    
    string ExpLog(long n)
    {
      return (1 + 1m / n).Pow(n).ToString();
    }
    
    void Out(string name, Func<long, string> func, long n)
    {
      var timer = Stopwatch.StartNew();
      var y = func(n);
      timer.Stop();
      Console.WriteLine("{0,-32} {1} {2}", y, timer.Elapsed, name);
    }
    
    void Run(int max)
    {
      for (var m = 0; m <= max; m++)
      {
        var n = (long)Math.Pow(10, m);
        Console.WriteLine(string.Format("- {0:D2}:{1:N0} ", m, n).PadRight(58, '-'));
        Out("Standard", Standard, n);
        Out("Direct", Direct, n);
        Out("Binary", Binary, n);
        Out("ExpLog", ExpLog, n);
      }
    }
  
    static void Main()
    {
      new Tester().Run(18);
    }
  }
}

這個程序使用四種方法來計算冪函數：

1. 第 10 至 16 行的 Standard 方法使用 BigInteger.Pow 方法來計算冪函數。這個計算結果（在有效數字范圍內）是准確值，作為其他方法的標准。
2. 第 18 至 24 行的 Direct 方法直接將 x 乘上 n 遍來計算冪函數，是最沒技術含量的暴力方法。時間復雜度是 O(N)。
3. 第 26 至 32 行的 Binary 方法將 n 視為二進制數，根據其為 1 的位來計算冪函數。這是經典的算法，時間復雜度是 O(logN)。FCL 的 BigInteger.Pow 方法也是使用這個算法。
4. 第 34 至 37 行的 ExpLog 方法使用 decimal 的擴展方法 Pow 來計算冪函數，是通過對數函數和指數函數來計算的：。理論上說，時間復雜度是 O(1)。

decimal 的擴展方法

下面就是 DecimalExtensions.cs：

using System;

namespace Skyiv.Extensions
{
  static class DecimalExtensions
  {
    static readonly int[] mask = { 1, 2, 4, 8, 16, 32, 64 };
    static readonly decimal ln10 = 2.3025850929940456840179914547m;
    static readonly decimal lnr = 0.2002433314278771112016301167m;
    static readonly decimal expmax = 66.542129333754749704054283659m;
    static readonly decimal[] exps =
    {
      2.71828182845904523536028747135m, // exp(1)
      7.38905609893065022723042746058m, // exp(2)
      54.5981500331442390781102612029m, // exp(4)
      2980.95798704172827474359209945m, // exp(8)
      8886110.52050787263676302374078m, // exp(16)
      78962960182680.6951609780226351m, // exp(32)
      6235149080811616882909238708.93m  // exp(64)
    };

    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.10527199m; k++) x /= 10;
      for (; x <= 0.1m; k--) x *= 10;        // ( 0.1000, 1.10527199 )
      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;
    }
    
    public static decimal Exp(this decimal x)
    {
      if (x > expmax) throw new OverflowException("overflow");
      if (x < -66) return 0;
      var n = (int)decimal.Round(x);
      if (n > 66) n--;
      decimal z = 1, y = Exponential(x - n);
      for (int m = (n < 0) ? -n : n, i = 0; i < mask.Length; i++)
        if ((m & mask[i]) != 0) z *= exps[i];
      return (n < 0) ? (y / z) : (y * z);
    }
    
    static decimal Exponential(decimal q)
    { // q (almost) in [ -0.5, 0.5 ]
      decimal y = 1, t = q;
      for (var i = 1; t != 0; t *= q / ++i) y += t;
      return y;
    }
    
    public static decimal Pow(this decimal x, decimal y)
    {
      if (x == 0 && y > 0) return 0;
      if (y == 0 && x != 0) return 1;
      return Exp(y * Log(x));
    }
  }
}

這個程序的詳細說明請見參考資料[5]和[6]。

編譯和運行

在 Arch Linux 操作系統的 Mono 環境下編譯和運行：

work$ dmcs -r:System.Numerics.dll Tester.cs DecimalExtensions.cs
work$ mono Tester.exe
- 00:1 ---------------------------------------------------
2.                               00:00:00.0085818 Standard
2                                00:00:00.0033230 Direct
2                                00:00:00.0002739 Binary
2.0000000000000000000000000005   00:00:00.0049157 ExpLog
- 01:10 --------------------------------------------------
2.5937424601                     00:00:00.0015421 Standard
2.5937424601000000000000000000   00:00:00.0000146 Direct
2.5937424601000000000000000000   00:00:00.0000092 Binary
2.5937424600999999999999999977   00:00:00.0000488 ExpLog
- 02:100 -------------------------------------------------
2.704813829421526093267194710807 00:00:00.0006872 Standard
2.7048138294215260932671947112   00:00:00.0000735 Direct
2.7048138294215260932671947103   00:00:00.0000234 Binary
2.7048138294215260932671947257   00:00:00.0000330 ExpLog
- 03:1,000 -----------------------------------------------
2.716923932235892457383088121947 00:00:00.0277308 Standard
2.7169239322358924573830881229   00:00:00.0007167 Direct
2.7169239322358924573830881218   00:00:00.0000159 Binary
2.7169239322358924573830883380   00:00:00.0000310 ExpLog
- 04:10,000 ----------------------------------------------
2.718145926825224864037664674913 00:00:03.3247007 Standard
2.7181459268252248640376646760   00:00:00.0068304 Direct
2.7181459268252248640376646665   00:00:00.0000191 Binary
2.7181459268252248640376679109   00:00:00.0000276 ExpLog
- 05:100,000 ---------------------------------------------
2.718268237174489668035064824426 00:07:56.2341075 Standard
2.7182682371744896680350648397   00:00:00.0686007 Direct
2.7182682371744896680350643783   00:00:00.0000222 Binary
2.7182682371744896680350286262   00:00:00.0000255 ExpLog
- 06:1,000,000 -------------------------------------------
Skip                             00:00:00.0000008 Standard
2.7182804693193768838197997202   00:00:00.6837104 Direct
2.7182804693193768838198166432   00:00:00.0000241 Binary
2.7182804693193768838199803836   00:00:00.0000213 ExpLog
- 07:10,000,000 ------------------------------------------
Skip                             00:00:00.0000009 Standard
2.7182816925449662711985502083   00:00:06.8334721 Direct
2.7182816925449662711985623547   00:00:00.0000289 Binary
2.7182816925449662712010419841   00:00:00.0000221 ExpLog
- 08:100,000,000 -----------------------------------------
Skip                             00:00:00.0000009 Standard
2.7182818148676362176529774118   00:01:08.3492423 Direct
2.7182818148676362176523859621   00:00:00.0000409 Binary
2.7182818148676362176710998015   00:00:00.0000230 ExpLog
- 09:1,000,000,000 ---------------------------------------
Skip                             00:00:00.0000007 Standard
2.7182818270999043223766453801   00:11:23.4187574 Direct
2.7182818270999043223770801045   00:00:00.0000442 Binary
2.7182818270999043220142064477   00:00:00.0000215 ExpLog
- 10:10,000,000,000 --------------------------------------
Skip                             00:00:00.0000007 Standard
Skip                             00:00:00.0000008 Direct
2.7182818283231311436196542093   00:00:00.0000349 Binary
2.7182818283231311439407330619   00:00:00.0000172 ExpLog
- 11:100,000,000,000 -------------------------------------
Skip                             00:00:00.0000008 Standard
Skip                             00:00:00.0000010 Direct
2.7182818284454538261539965115   00:00:00.0000398 Binary
2.7182818284454538262180262237   00:00:00.0000176 ExpLog
- 12:1,000,000,000,000 -----------------------------------
Skip                             00:00:00.0000010 Standard
Skip                             00:00:00.0000007 Direct
2.7182818284576860942863185484   00:00:00.0000403 Binary
2.7182818284576860944460582886   00:00:00.0000174 ExpLog
- 13:10,000,000,000,000 ----------------------------------
Skip                             00:00:00.0000009 Standard
Skip                             00:00:00.0000007 Direct
2.7182818284589093212295138270   00:00:00.0000436 Binary
2.7182818284589093212688645227   00:00:00.0000176 ExpLog
- 14:100,000,000,000,000 ---------------------------------
Skip                             00:00:00.0000009 Standard
Skip                             00:00:00.0000009 Direct
2.7182818284590316438350187680   00:00:00.0000480 Binary
2.7182818284590452353602874714   00:00:00.0000112 ExpLog
- 15:1,000,000,000,000,000 -------------------------------
Skip                             00:00:00.0000009 Standard
Skip                             00:00:00.0000009 Direct
2.7182818284590431765145511000   00:00:00.0000522 Binary
2.7182818284590452353602874714   00:00:00.0000114 ExpLog
- 16:10,000,000,000,000,000 ------------------------------
Skip                             00:00:00.0000009 Standard
Skip                             00:00:00.0000006 Direct
2.7182818284590335325626228124   00:00:00.0000547 Binary
2.7182818284590452353602874714   00:00:00.0000109 ExpLog
- 17:100,000,000,000,000,000 -----------------------------
Skip                             00:00:00.0000010 Standard
Skip                             00:00:00.0000006 Direct
2.7182818284590296936415060358   00:00:00.0000567 Binary
2.7182818284590452353602874714   00:00:00.0000108 ExpLog
- 18:1,000,000,000,000,000,000 ---------------------------
Skip                             00:00:00.0000010 Standard
Skip                             00:00:00.0000006 Direct
2.7182818284590434884535909399   00:00:00.0000615 Binary
2.7182818284590452353602874714   00:00:00.0000108 ExpLog
work$ echo 'scale=30;e(1)' | bc -lq
2.718281828459045235360287471352

在上述結果中：

1. 最后一行是 e 的近似值，是使用 Linux 操作系統的高精度計算器 bc 計算的，請見參考資料[4]。
2. 使用 BigInteger.Pow 計算出來的是准確值。在 05:100,000 這一組中，計算結果達 500,001 個十進制數字。當 n 達到 10⁶ 以后，由於計算量太大，已經無法在合理的時間內計算准確值了。
3. 使用 Direct 計算最慢（除 Standard 外，因為計算量不同）。當 n 達到 10¹⁰ 以后，由於費時太多，已經不使用 Direct 方法計算了。
4. 使用 Binary 計算的速度非常快，其精度和 Direct 差不多。這兩者答案不同說明 decimal 的乘法不滿足結合律。
5. 使用 ExpLog 計算的速度理論上是最快的，實際的速度和 ExpLog 差不多，因為 n 還不夠大。其精度在 n 不是很大時稍差。
6. 當 n 達到 10¹⁴ 以后，ExpLog 計算出來的值在 29 個有效數字范圍內已經等於 e 值，不再變化了。
7. 當 n 達到 10¹⁴ 以后，由於舍入誤差的累計，Binary 計算出來的值大約只有 14 個有效數字是可信的，再增大 n 值也不能更逼近 e 值了。也就是說，在逼近 e 值的意義上說，計算結果在有效數字范圍內不再變化了。
8. 要計算的冪函數是增函數，請注意觀察上述運行結果是如何體現這一點的。

參數資料

1. Wikipedia: e (mathematical constant)
2. MSDN: BigInteger.Pow 方法 (System.Numerics)
3. MSDN: Math.Pow 方法 (System)
4. Linux man pages: bc - An arbitrary precision calculator language
5. 博客園: 計算自然對數的算法
6. 博客園: 計算指數函數的算法