计算机指数函数算法框图,计算指数函数的算法

引言

我在上一篇随笔中介绍了计算自然对数的快速算法。现在我们来看看计算指数函数的算法。我们知道&＃xff0c;指数函数 ex 可以展开为泰勒级数&＃xff1a;

这个级数对全体实数 x 都收敛&＃xff0c;并且在 x 接近零时收敛得比较快。

实现该算法的 C# 程序

根据前面所述的 ex 的泰勒级数展开式&＃xff0c;可以写出以下 C# 程序来为 decimal 数据类型添加一个 Exp 扩展方法&＃xff1a;

1 usingSystem;2

3 namespaceSkyiv.Extensions4 {5 static classDecimalExtensions6 {7 static readonly decimal expmax &＃61; 66.542129333754749704054283659m;8 static readonly int[] mask &＃61; { 1, 2, 4, 8, 16, 32, 64};9 static readonly decimal[] exps &＃61;

10 {11 2.71828182845904523536028747135m, //exp(1)

12 7.38905609893065022723042746058m, //exp(2)

13 54.5981500331442390781102612029m, //exp(4)

14 2980.95798704172827474359209945m, //exp(8)

15 8886110.52050787263676302374078m, //exp(16)

16 78962960182680.6951609780226351m, //exp(32)

17 6235149080811616882909238708.93m //exp(64)

18 };19

20 public static decimal Exp(this decimalx)21 {22 if (x > expmax) throw new OverflowException("overflow");23 if (x <-66) return 0;24 var n &＃61; (int)decimal.Round(x);25 if (n > 66) n--;26 decimal z &＃61; 1, y &＃61; Exponential(x -n);27 for (int m &＃61; (n <0) ? -n : n, i &＃61; 0; i

32 static decimal Exponential(decimalq)33 { //q (almost) in [ -0.5, 0.5 ]

34 decimal y &＃61; 1, t &＃61;q;35 for (var i &＃61; 1; t !&＃61; 0; t *&＃61; q / &＃43;&＃43;i) y &＃43;&＃61;t;36 returny;37 }38 }39 }

简要说明如下&＃xff1a;

第 7 行的 expmax 的值是 decimal.MaxValue 的自然对数的近似值&＃xff0c;用于检测 Exp 方法是否溢出(第 22 行)。

第 20 至 30 行的 Exp 扩展方法就是用来计算指数函数了。

该方法利用 ex&＃43;y &＃61; exey 这个公式&＃xff0c;将参数 x 分为整数部分 n 和小数部分 x - n 来计算。

整数部分 n 又分解为 1、2、4、8、16、32、 64 诸数中某些的和&＃xff0c;利用事先计算出来的常量来计算。

第 25 行是为了防止将 e66.5421 分解为 e67e-0.4579&＃xff0c;这样在计算 e67 时会溢出。而是分解为 e66e0.5421。

第 32 至 37 行的 Exponential 方法使用泰勒级数来计算 ex 。它的参数 q 越接近于零就计算得越快。

这个算法还是很快的&＃xff0c;第 35 行的 for 循环执行次数不会超过 22 次。

测试程序

下面就是调用 decimal 数据类型的 Exp 扩展方法的测试程序&＃xff1a;

1 usingSystem;2 usingSkyiv.Extensions;3

4 classTester5 {6 static voidMain()7 {8 try

9 {10 foreach (var x in new decimal[] {11 -100, -66, -65, -1, 0, 1, 2.5m, 16, 66.5421m, 67})12 Console.WriteLine("{0,-30}: exp({1})", x.Exp(), x);13 }14 catch(Exception ex) { Console.WriteLine(ex.Message); }15 }16 }

运行结果如下所示&＃xff1a;

work$ dmcs Tester.cs DecimalExtensions.cs

work$ mono Tester.exe

0 : exp(-100)

0.0000000000000000000000000000: exp(-66)

0.0000000000000000000000000001: exp(-65)

0.3678794411714423215955237702: exp(-1)

1 : exp(0)

2.7182818284590452353602874714: exp(1)

12.182493960703473438070175950: exp(2.5)

8886110.520507872636763023741 : exp(16)

79225838488862236701995526357 : exp(66.5421)

overflow

可以看出&＃xff0c;在计算 e67 时发现了溢出。这是因为&＃xff1a;

decimal.MaxValue &＃61; 79,228,162,514,264,337,593,543,950,335

e67 &＃61; 125,236,317,084,221,378,051,352,196,074.4365767534885274 ...

可以看出&＃xff0c;e67 已经超过 decimal 的最大值了。

事先计算的常数

在 DecimalExtensions.cs 程序的第 9 至 18 行中的 exps 静态只读数组中存放了 e1、e2、e4、e8、e16、e32 和 e64 的值。这些值是如何得到的呢&＃xff1f;这很简单&＃xff0c;Linux 操作系统中有一个高精度计算器 bc 。我们可以先编辑一个如下内容的文本文件 exps_in.txt&＃xff1a;

scale&＃61;30

e(1)

e(2)

e(4)

e(8)

e(16)

e(32)

e(64)

l(2^96-1)

quit

上面的 e 代表 exp&＃xff0c;l 代表 ln&＃xff0c;296 - 1 就是 decimal.MaxValue。然后执行以下命令&＃xff1a;

work$ bc -l exps_in.txt > exps_out.txt

就可以得出如下内容的输出 exps_out.txt&＃xff1a;

2.718281828459045235360287471352

7.389056098930650227230427460575

54.598150033144239078110261202860

2980.957987041728274743592099452888

8886110.520507872636763023740781450350

78962960182680.695160978022635108224219956195

6235149080811616882909238708.928469744831391846235799914388

66.542129333754749704054283659972

稍加整理&＃xff0c;就可以用在上述 C# 程序中了&＃xff1a;

前 7 行就是 e 的各次幂。

最后一行就是 decimal.MaxValue 的自然对数。

参考资料