计算幂函数的几种方法
来源:互联网 发布:淘宝一键生成手机详情 编辑:程序博客网 时间:2024/05/17 02:09
引言
我们知道,自然对数的底 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^m12 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 = Stopwatch.StartNew();42 var y = func(n);43 timer.Stop();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.cswork$ mono Tester.exe- 00:1 ---------------------------------------------------2. 00:00:00.0085818 Standard2 00:00:00.0033230 Direct2 00:00:00.0002739 Binary2.0000000000000000000000000005 00:00:00.0049157 ExpLog- 01:10 --------------------------------------------------2.5937424601 00:00:00.0015421 Standard2.5937424601000000000000000000 00:00:00.0000146 Direct2.5937424601000000000000000000 00:00:00.0000092 Binary2.5937424600999999999999999977 00:00:00.0000488 ExpLog- 02:100 -------------------------------------------------2.704813829421526093267194710807 00:00:00.0006872 Standard2.7048138294215260932671947112 00:00:00.0000735 Direct2.7048138294215260932671947103 00:00:00.0000234 Binary2.7048138294215260932671947257 00:00:00.0000330 ExpLog- 03:1,000 -----------------------------------------------2.716923932235892457383088121947 00:00:00.0277308 Standard2.7169239322358924573830881229 00:00:00.0007167 Direct2.7169239322358924573830881218 00:00:00.0000159 Binary2.7169239322358924573830883380 00:00:00.0000310 ExpLog- 04:10,000 ----------------------------------------------2.718145926825224864037664674913 00:00:03.3247007 Standard2.7181459268252248640376646760 00:00:00.0068304 Direct2.7181459268252248640376646665 00:00:00.0000191 Binary2.7181459268252248640376679109 00:00:00.0000276 ExpLog- 05:100,000 ---------------------------------------------2.718268237174489668035064824426 00:07:56.2341075 Standard2.7182682371744896680350648397 00:00:00.0686007 Direct2.7182682371744896680350643783 00:00:00.0000222 Binary2.7182682371744896680350286262 00:00:00.0000255 ExpLog- 06:1,000,000 -------------------------------------------Skip 00:00:00.0000008 Standard2.7182804693193768838197997202 00:00:00.6837104 Direct2.7182804693193768838198166432 00:00:00.0000241 Binary2.7182804693193768838199803836 00:00:00.0000213 ExpLog- 07:10,000,000 ------------------------------------------Skip 00:00:00.0000009 Standard2.7182816925449662711985502083 00:00:06.8334721 Direct2.7182816925449662711985623547 00:00:00.0000289 Binary2.7182816925449662712010419841 00:00:00.0000221 ExpLog- 08:100,000,000 -----------------------------------------Skip 00:00:00.0000009 Standard2.7182818148676362176529774118 00:01:08.3492423 Direct2.7182818148676362176523859621 00:00:00.0000409 Binary2.7182818148676362176710998015 00:00:00.0000230 ExpLog- 09:1,000,000,000 ---------------------------------------Skip 00:00:00.0000007 Standard2.7182818270999043223766453801 00:11:23.4187574 Direct2.7182818270999043223770801045 00:00:00.0000442 Binary2.7182818270999043220142064477 00:00:00.0000215 ExpLog- 10:10,000,000,000 --------------------------------------Skip 00:00:00.0000007 StandardSkip 00:00:00.0000008 Direct2.7182818283231311436196542093 00:00:00.0000349 Binary2.7182818283231311439407330619 00:00:00.0000172 ExpLog- 11:100,000,000,000 -------------------------------------Skip 00:00:00.0000008 StandardSkip 00:00:00.0000010 Direct2.7182818284454538261539965115 00:00:00.0000398 Binary2.7182818284454538262180262237 00:00:00.0000176 ExpLog- 12:1,000,000,000,000 -----------------------------------Skip 00:00:00.0000010 StandardSkip 00:00:00.0000007 Direct2.7182818284576860942863185484 00:00:00.0000403 Binary2.7182818284576860944460582886 00:00:00.0000174 ExpLog- 13:10,000,000,000,000 ----------------------------------Skip 00:00:00.0000009 StandardSkip 00:00:00.0000007 Direct2.7182818284589093212295138270 00:00:00.0000436 Binary2.7182818284589093212688645227 00:00:00.0000176 ExpLog- 14:100,000,000,000,000 ---------------------------------Skip 00:00:00.0000009 StandardSkip 00:00:00.0000009 Direct2.7182818284590316438350187680 00:00:00.0000480 Binary2.7182818284590452353602874714 00:00:00.0000112 ExpLog- 15:1,000,000,000,000,000 -------------------------------Skip 00:00:00.0000009 StandardSkip 00:00:00.0000009 Direct2.7182818284590431765145511000 00:00:00.0000522 Binary2.7182818284590452353602874714 00:00:00.0000114 ExpLog- 16:10,000,000,000,000,000 ------------------------------Skip 00:00:00.0000009 StandardSkip 00:00:00.0000006 Direct2.7182818284590335325626228124 00:00:00.0000547 Binary2.7182818284590452353602874714 00:00:00.0000109 ExpLog- 17:100,000,000,000,000,000 -----------------------------Skip 00:00:00.0000010 StandardSkip 00:00:00.0000006 Direct2.7182818284590296936415060358 00:00:00.0000567 Binary2.7182818284590452353602874714 00:00:00.0000108 ExpLog- 18:1,000,000,000,000,000,000 ---------------------------Skip 00:00:00.0000010 StandardSkip 00:00:00.0000006 Direct2.7182818284590434884535909399 00:00:00.0000615 Binary2.7182818284590452353602874714 00:00:00.0000108 ExpLogwork$ echo 'scale=30;e(1)' | bc -lq2.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 值的意义上说,计算结果在有效数字范围内不再变化了。
- 要计算的幂函数是增函数,请注意观察上述运行结果是如何体现这一点的。
0 0
- 计算幂函数的几种方法
- 计算幂函数的几种方法
- 几种计算时间的方法
- 计算fibnacci 级数的几种方法
- 计算组合数的几种方法
- 计算星期几的方法
- 计算一个数字的长度的几种方法
- 计算两个数的平均数的几种方法解读
- javascript调用函数的几种方法
- Matlab自定义函数的几种方法
- strlen函数实现的几种方法
- Matlab自定义函数的几种方法
- scala函数定义的几种方法
- Matlab自定义函数的几种方法
- Matlab自定义函数的几种方法
- js定义函数的几种方法
- strlen函数的几种实现方法
- 求欧拉函数的几种方法
- 【Qt OpenGL教程】01:创建一个OpenGL窗口
- c# winform DataGridView选择一整行的相关属性
- Sql常见函数大全
- SQL语句笔记
- uva 11646
- 计算幂函数的几种方法
- SQLServer权限
- poj_2478_欧拉函数
- html/css的网页制作
- 两种计算自然对数的算法比较
- 计算指数函数的算法
- C++ 管道
- 浮点数内存表示方法
- Activity数据传输--Intent显示、隐式、不同应用之间启动