计算指数函数的算法

Posted zhchoutai

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了计算指数函数的算法相关的知识,希望对你有一定的参考价值。

引言

我在上一篇随笔中介绍了计算自然对数的高速算法。如今我们来看看计算指数函数的算法。我们知道。指数函数 ex 能够展开为泰勒级数:

exp

这个级数对全体实数 x 都收敛,而且在 x 接近零时收敛得比較快。

实现该算法的 C# 程序

依据前面所述的 ex 的泰勒级数展开式,能够写出下面 C# 程序来为 decimal 数据类型加入一个 Exp 扩展方法:

复制代码
 1 using System;
 2 
 3 namespace Skyiv.Extensions
 4 {
 5   static class DecimalExtensions
 6   {
 7     static readonly decimal expmax = 66.542129333754749704054283659m;
 8     static readonly int[] mask = { 1, 2, 4, 8, 16, 32, 64 };
 9     static readonly decimal[] exps =
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 decimal x)
21     {
22       if (x > expmax) throw new OverflowException("overflow");
23       if (x < -66) return 0;
24       var n = (int)decimal.Round(x);
25       if (n > 66) n--;
26       decimal z = 1, y = Exponential(x - n);
27       for (int m = (n < 0) ?

-n : n, i = 0; i < mask.Length; i++) 28 if ((m & mask[i]) != 0) z *= exps[i]; 29 return (n < 0) ? (y / z) : (y * z); 30 } 31 32 static decimal Exponential(decimal q) 33 { // q (almost) in [ -0.5, 0.5 ] 34 decimal y = 1, t = q; 35 for (var i = 1; t != 0; t *= q / ++i) y += t; 36 return y; 37 } 38 } 39 }

复制代码

简要说明例如以下:

  1. 第 7 行的 expmax 的值是 decimal.MaxValue 的自然对数的近似值,用于检測 Exp 方法是否溢出(第 22 行)。

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

  3. 该方法利用 ex+y = exey 这个公式。将參数 x 分为整数部分 n 和小数部分 x - n 来计算。
  4. 整数部分 n 又分解为 1、2、4、8、16、32、 64 诸数中某些的和,利用事先计算出来的常量来计算。
  5. 第 25 行是为了防止将 e66.5421 分解为 e67e-0.4579。这样在计算 e67 时会溢出。而是分解为 e66e0.5421
  6. 第 32 至 37 行的 Exponential 方法使用泰勒级数来计算 ex 。它的參数 q 越接近于零就计算得越快。

  7. 这个算法还是非常快的,第 35 行的 for 循环运行次数不会超过 22 次。

測试程序

以下就是调用 decimal 数据类型的 Exp 扩展方法的測试程序:

复制代码
 1 using System;
 2 using Skyiv.Extensions;
 3 
 4 class Tester
 5 {
 6   static void Main()
 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 }
复制代码

执行结果例如以下所看到的:

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

能够看出,在计算 e67 时发现了溢出。这是由于:

  • decimal.MaxValue = 79,228,162,514,264,337,593,543,950,335
  • e67 = 125,236,317,084,221,378,051,352,196,074.4365767534885274 ...

能够看出,e67 已经超过 decimal 的最大值了。

事先计算的常数

在 DecimalExtensions.cs 程序的第 9 至 18 行中的 exps 静态仅仅读数组中存放了 e1、e2、e4、e8、e16、e32 和 e64 的值。这些值是怎样得到的呢?这非常easy,Linux 操作系统中有一个高精度计算器 bc 。

我们能够先编辑一个例如以下内容的文本文件 exps_in.txt:

scale=30
e(1)
e(2)
e(4)
e(8)
e(16)
e(32)
e(64)
l(2^96-1)
quit

上面的 e 代表 exp。l 代表 ln,296 - 1 就是 decimal.MaxValue。

然后运行下面命令:

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

就能够得出例如以下内容的输出 exps_out.txt:

2.718281828459045235360287471352
7.389056098930650227230427460575
54.598150033144239078110261202860
2980.957987041728274743592099452888
8886110.520507872636763023740781450350
78962960182680.695160978022635108224219956195
6235149080811616882909238708.928469744831391846235799914388
66.542129333754749704054283659972

稍加整理,就能够用在上述 C# 程序中了:

  • 前 7 行就是 e 的各次幂。

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

以上是关于计算指数函数的算法的主要内容,如果未能解决你的问题,请参考以下文章

视频二次裁剪时间计算出原片的时间片段算法

10个JavaScript代码片段,使你更加容易前端开发。

10个JavaScript代码片段,使你更加容易前端开发。

以下代码片段的算法复杂度

计算机视觉常用算法讲解

有人可以解释啥是 SVN 平分算法吗?理论上和通过代码片段[重复]