计算幂函数的几种方法

引言

我们知道，自然对数的底 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 = 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 }

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

  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：

 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

在上述结果中：

  1. 最后一行是 e 的近似值，是使用 Linux 操作系统的高精度计算器 bc 计算的，请见参考资料[4]。
  2. 使用 BigInteger.Pow 计算出来的是准确值。在 05:100,000 这一组中，计算结果达 500,001 个十进制数字。当 n 达到 106 以后，由于计算量太大，已经无法在合理的时间内计算准确值了。
  3. 使用 Direct 计算最慢（除 Standard 外，因为计算量不同）。当 n 达到 1010 以后，由于费时太多，已经不使用 Direct 方法计算了。
  4. 使用 Binary 计算的速度非常快，其精度和 Direct 差不多。这两者答案不同说明 decimal 的乘法不满足结合律。
  5. 使用 ExpLog 计算的速度理论上是最快的，实际的速度和 ExpLog 差不多，因为 n 还不够大。其精度在 n 不是很大时稍差。
  6. 当 n 达到 1014 以后，ExpLog 计算出来的值在 29 个有效数字范围内已经等于 e 值，不再变化了。
  7. 当 n 达到 1014 以后，由于舍入误差的累计，Binary 计算出来的值大约只有 14 个有效数字是可信的，再增大 n 值也不能更逼近 e 值了。也就是说，在逼近 e 值的意义上说，计算结果在有效数字范围内不再变化了。
  8. 要计算的幂函数是增函数，请注意观察上述运行结果是如何体现这一点的。

版权声明：本文为博主http://www.zuiniusn.com原创文章，未经博主允许不得转载。

时间： 2024-11-15 08:15:19

计算幂函数的几种方法的相关文章

获取时间和计算时间差的几种方法总结

转载自:http://blog.csdn.net/coder_xia/article/details/6566708 一.标准C和C++都可用 1.获取时间用time_t time( time_t * timer ),计算时间差使用double difftime( time_t timer1, time_t timer0 ). 精确到秒. 测试程序如下: 1 #include <time.h> 2 #include <stdio.h> 3 int main() 4 { 5 time

【转载】c/c++在windows下获取时间和计算时间差的几种方法总结

一.标准C和C++都可用 1.获取时间用time_t time( time_t * timer ),计算时间差使用double difftime( time_t timer1, time_t timer0 ). 精确到秒. 测试程序如下: #include <time.h> #include <stdio.h> int main() { time_t start ,end ; double cost; time(&start); sleep(1); time(&en

C#中Math类的计算整数的三种方法

1.Math.Round:四舍六入五取偶 引用内容 Math.Round(0.0) //0 Math.Round(0.1) //0 Math.Round(0.2) //0 Math.Round(0.3) //0 Math.Round(0.4) //0 Math.Round(0.5) //0 Math.Round(0.6) //1 Math.Round(0.7) //1 Math.Round(0.8) //1 Math.Round(0.9) //1 说明:对于1.5,因要返回偶数,所以结果为2.

康复计划#3 简单常用的几种计算自然数幂和的方法

本篇口胡写给我自己这样的东西都忘光的残废选手 以及暂时还不会自然数幂和的人- 这里大概给出最简单的几种方法:扰动法(化为递推式),斯特林数(离散微积分),高阶差分(牛顿级数),伯努利数(指数生成函数)- 不同方法的思维难度.普适程度.实现难度.时间复杂度上面都有差异-同时自然数幂和是探究各种求和方法的经典例子,了解多一点它的做法对于处理各种求和问题是有所帮助的- 问题:求$\sum_{k=0}^{n} k^t$,其中$t \in \mathbb{N}$是一个常数.要求求解的时间复杂度与$n$无关

Java 计算中英文长度的若干种方法

在项目开发中经常碰到到输入字符的校验,特别是中英文混合在一起的校验.而为了满足校验的需求,有时需要计算出中英文的长度. 本文将通过几种常用的方法实现长度的计算: <span style="font-size:18px;">import java.io.UnsupportedEncodingException; /** * 中英文校验的处理 * @author a123demi * */ public class EnChValidate { public static vo

【九天教您南方cass 9.1】 10 DTM土方计算的四种方法

同学们大家好,欢迎收看由老王测量上班记出品的cass9.1视频课程 我是本节课主讲老师九天. 我们讲课的教程附件也是共享的,请注意索取测量空间中. 九天老师的联系方式  点击直接请教九天老师吧! DTM土方计算的四种方法 课程教案 一.根据坐标文件. (1).-->根据坐标文件 (2).-->指定坐标文件 (3).-->DTM土方计算参数设置 (4).-->得出结果dtmtf.log (5).-->并指定三角网法土石方计算成果表的位置 二.根据图上高程点. (1).-->

聊聊JVM（三）两种计算Java对象大小的方法

这篇说说如何计算Java对象大小的方法.之前在聊聊高并发(四)Java对象的表示模型和运行时内存表示 这篇中已经说了Java对象的内存表示模型是Oop-Klass模型. 普通对象的结构如下,按64位机器的长度计算 1. 对象头(_mark), 8个字节 2. Oop指针,如果是32G内存以下的,默认开启对象指针压缩,4个字节 3. 数据区 4.Padding(内存对齐),按照8的倍数对齐 数组对象结构是 1. 对象头(_mark), 8个字节 2. Oop指针,如果是32G内存以下的,默认开启对

Linux下计算字符串长度的四种方法

在linux运维场景中,经常会碰到计算字符长度的场景,这里罗列四种方法: [[email protected] ~]# char="I love you" 方法一:[[email protected] ~]# echo ${#char} 方法二:[[email protected] ~]# expr length "$char" 方法三:[[email protected] ~]# echo $char|wc -L 方法四:[[email protected] ~]

26 计算用户输入的内容中索引为奇数并且对应的元素为数字的个数的两种方法

#计算用户输入的内容中索引为奇数并且对应的元素为数字的个数第二种方法content = input(">>>")count = 0for i in range(len(content)):#i就是下标,或者说就是索引 if i % 2 == 1 and content[i].isdigit(): count += 1print(count) #计算用户输入的内容中索引为奇数并且对应的元素为数字的个数的第一种方法 li = []res = ' '.join(input(