计算幂函数的几种方法 武汉北大青鸟技术介绍

 武汉北大青鸟技术介绍计算幂函数的几种方法。

    我们知道，自然对数的底 e 定义为以下限值：

　　这个公式很适合于对幂函数的计算进行一些测试，得到的结果是 e 的近似值，不用担心当 n 很大时计算结果会溢出。

　　测试程序

　　下面就是 Tester.cs：

　　1 using System;

　　2 using System.Numerics;

　　3 using System.Diagnostics;

　　4 using Skyiv.Extensions;

　　5

　　6 namespace Skyiv.Test

　　7 {

　　8   sealed class Tester

　　9   {

　　10     string Standard(long n)

　　11     { // n == 10^m

　　12       if (n > 100000) return "Skip";

　　13       var s = BigInteger.Pow(n + 1, (int)n).ToString();

　　14       s = s.Substring(0, Math.Min(31, s.Length));

　　15       return s[0] + "." + s.Substring(1);

　　16     }

　　17

　　18     string Direct(long n)

　　19     {

　　20       if (n > 1000000000) return "Skip";

　　21       var y = 1m;

　　22       for (var x = 1 + 1m / n; n > 0; n--) y *= x;

　　23       return y.ToString();

　　24     }

　　25

　　26     string Binary(long n)

　　27     {

　　28       var y = 1m;

　　29       for (var x = 1 + 1m / n; n != 0; x *= x, n >>= 1)

　　30         if ((n & 1) != 0) y *= x;

　　31       return y.ToString();

　　32     }

　　33

　　34     string ExpLog(long n)

　　35     {

　　36       return (1 + 1m / n).Pow(n).ToString();

　　37     }

　　38

　　39     void Out(string name, Func<long, string> func, long n)

　　40     {

　　41       var timer = Swatch.StartNew();

　　42       var y = func(n);

　　43       timer.S();

　　44       Console.WriteLine("{0,-32} {1} {2}", y, timer.Elapsed, name);

　　45     }

　　46

　　47     void Run(int max)

　　48     {

　　49       for (var m = 0; m <= max; m++)

　　50       {

　　51         var n = (long)Math.Pow(10, m);

　　52         Console.WriteLine(string.Format("- {0:D2}:{1:N0} ", m, n).PadRight(58, '-'));

　　53         Out("Standard", Standard, n);

　　54         Out("Direct", Direct, n);

　　55         Out("Binary", Binary, n);

　　56         Out("ExpLog", ExpLog, n);

　　57       }

　　58     }

　　59

　　60     static void Main()

　　61     {

　　62       new Tester().Run(18);

　　63     }

　　64   }

　　65 }

　这个程序使用四种方法来计算幂函数：

　　第 10 至 16 行的 Standard 方法使用 BigInteger.Pow 方法来计算幂函数。这个计算结果（在有效数字范围内）是准确值，作为其他方法的标准。

　　第 18 至 24 行的 Direct 方法直接将 x 乘上 n 遍来计算幂函数，是没技术含量的暴力方法。时间复杂度是 O(N)。

　　第 26 至 32 行的 Binary 方法将 n 视为二进制数，根据其为 1 的位来计算幂函数。这是经典的算法，时间复杂度是 O(logN)。FCL 的 BigInteger.Pow 方法也是使用这个算法。

　　第 34 至 37 行的 ExpLog 方法使用 decimal 的扩展方法 Pow 来计算幂函数，是通过对数函数和指数函数来计算的：。理论上说，时间复杂度是 O(1)。

　　decimal 的扩展方法

　　下面就是 DecimalExtensions.cs：

　　1 using System;

　　2

　　3 namespace Skyiv.Extensions

　　4 {

　　5   static class DecimalExtensions

　　6   {

　　7     static readonly int[] mask = { 1, 2, 4, 8, 16, 32, 64 };

　　8     static readonly decimal ln10 = 2.3025850929940456840179914547m;

　　9     static readonly decimal lnr = 0.2002433314278771112016301167m;

　　10     static readonly decimal expmax = 66.542129333754749704054283659m;

　　11     static readonly decimal[] exps =

　　12     {

　　13       2.71828182845904523536028747135m, // exp(1)

　　14       7.38905609893065022723042746058m, // exp(2)

　　15       54.5981500331442390781102612029m, // exp(4)

　　16       2980.95798704172827474359209945m, // exp(8)

　　17       8886110.52050787263676302374078m, // exp(16)

　　18       78962960182680.6951609780226351m, // exp(32)

　　19       6235149080811616882909238708.93m  // exp(64)

　　20     };

　　21

　　22     public static decimal Log10(this decimal x)

　　23     {

　　24       return Log(x) / ln10;

　　25     }

　　26

　　27     public static decimal Log(this decimal x)

　　28     {

　　29       if (x <= 0) throw new ArgumentException("Must be positive");

　　30       int k = 0, l = 0;

　　31       for (; x >= 1.10527199m; k++) x /= 10;

　　32       for (; x <= 0.1m; k--) x *= 10;        // ( 0.1000, 1.10527199 )

　　33       for (; x < 0.9047m; l--) x *= 1.2217m; // [ 0.9047, 1.10527199 )

　　34       return k * ln10 + l * lnr + Logarithm((x - 1) / (x + 1));

　　35     }

　　36

　　37     static decimal Logarithm(decimal y)

　　38     { // y in ( -0.05-, 0.05+ ), return ln((1+y)/(1-y))

　　39       decimal v = 1, y2 = y * y, t = y2, z = t / 3;

　　40       for (var i = 3; z != 0; z = (t *= y2) / (i += 2)) v += z;

　　41       return v * y * 2;

　　42     }

　　43

　　44     public static decimal Exp(this decimal x)

　　45     {

　　46       if (x > expmax) throw new OverflowException("overflow");

　　47       if (x < -66) return 0;

　　48       var n = (int)decimal.Round(x);

　　49       if (n > 66) n--;

　　50       decimal z = 1, y = Exponential(x - n);

　　51       for (int m = (n < 0) ? -n : n, i = 0; i < mask.Length; i++)

　　52         if ((m & mask[i]) != 0) z *= exps[i];

　　53       return (n < 0) ? (y / z) : (y * z);

　　54     }

　　55

　　56     static decimal Exponential(decimal q)

　　57     { // q (almost) in [ -0.5, 0.5 ]

　　58       decimal y = 1, t = q;

　　59       for (var i = 1; t != 0; t *= q / ++i) y += t;

　　60       return y;

　　61     }

　　62

　　63     public static decimal Pow(this decimal x, decimal y)

　　64     {

　　65       if (x == 0 && y > 0) return 0;

　　66       if (y == 0 && x != 0) return 1;

　　67       return Exp(y * Log(x));

　　68     }

　　69   }

　　70 }

　　这个程序的详细说明请见参考资料[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

　　在上述结果中：

　　后一行是 e 的近似值，是使用 Linux 操作系统的高精度计算器 bc 计算的，请见参考资料[4]。

　　使用 BigInteger.Pow 计算出来的是准确值。在 05:100,000 这一组中，计算结果达 500,001 个十进制数字。当 n 达到 106 以后，由于计算量太大，已经无法在合理的时间内计算准确值了。

　　使用 Direct 计算慢（除 Standard 外，因为计算量不同）。当 n 达到 1010 以后，由于费时太多，已经不使用 Direct 方法计算了。

　　使用 Binary 计算的速度非常快，其精度和 Direct 差不多。这两者答案不同说明 decimal 的乘法不满足结合律。

　　使用 ExpLog 计算的速度理论上是快的，实际的速度和 ExpLog 差不多，因为 n 还不够大。其精度在 n 不是很大时稍差。

　　当 n 达到 1014 以后，ExpLog 计算出来的值在 29 个有效数字范围内已经等于 e 值，不再变化了。

　　当 n 达到 1014 以后，由于舍入误差的累计，Binary 计算出来的值大约只有 14 个有效数字是可信的，再增大 n 值也不能更逼近 e 值了。也就是说，在逼近 e 值的意义上说，计算结果在有效数字范围内不再变化了。

　　要计算的幂函数是增函数，请注意观察上述运行结果是如何体现这一点的。

        了解北大青鸟武汉宏鹏光谷校区新技术请继续关注武汉北大青鸟官网。