计算指数函数的算法
引言
我在上一篇随笔中介绍了计算自然对数的快速算法。现在我们来看看计算指数函数的算法。我们知道,指数函数 ex可以展开为泰勒级数:
这个级数对全体实数 x 都收敛,并且在 x 接近零时收敛得比较快。
实现该算法的 C# 程序
根据前面所述的 ex的泰勒级数展开式,可以写出以下 C# 程序来为 decimal 数据类型添加一个 Exp 扩展方法:
using System;
namespace Skyiv.Extensions
{
static class DecimalExtensions
{
static readonly decimal expmax = 66.542129333754749704054283659m;
static readonly int[] mask = { 1, 2, 4, 8, 16, 32, 64 };
static readonly decimal[] exps =
{
2.71828182845904523536028747135m, // exp(1)
7.38905609893065022723042746058m, // exp(2)
54.5981500331442390781102612029m, // exp(4)
2980.95798704172827474359209945m, // exp(8)
8886110.52050787263676302374078m, // exp(16)
78962960182680.6951609780226351m, // exp(32)
6235149080811616882909238708.93m // exp(64)
};
public static decimal Exp(this decimal x)
{
if (x > expmax) throw new OverflowException("overflow");
if (x < -66) return 0;
var n = (int)decimal.Round(x);
if (n > 66) n--;
decimal z = 1, y = Exponential(x - n);
for (int m = (n < 0) ? -n : n, i = 0; i < mask.Length; i++)
if ((m & mask[i]) != 0) z *= exps[i];
return (n < 0) ? (y / z) : (y * z);
}
static decimal Exponential(decimal q)
{ // q (almost) in [ -0.5, 0.5 ]
decimal y = 1, t = q;
for (var i = 1; t != 0; t *= q / ++i) y += t;
return y;
}
}
}
简要说明如下:
- 第 7 行的 expmax 的值是 decimal.MaxValue 的自然对数的近似值,用于检测 Exp 方法是否溢出(第 22 行)。
- 第 20 至 30 行的 Exp 扩展方法就是用来计算指数函数了。
- 该方法利用 ex+y= exey这个公式,将参数 x 分为整数部分 n 和小数部分 x - n 来计算。
- 整数部分 n 又分解为 1、2、4、8、16、32、 64 诸数中某些的和,利用事先计算出来的常量来计算。
- 第 25 行是为了防止将 e66.5421分解为 e67e-0.4579,这样在计算 e67时会溢出。而是分解为 e66e0.5421。
- 第 32 至 37 行的 Exponential 方法使用泰勒级数来计算 ex。它的参数 q 越接近于零就计算得越快。
- 这个算法还是很快的,第 35 行的 for 循环执行次数不会超过 22 次。
测试程序
下面就是调用 decimal 数据类型的 Exp 扩展方法的测试程序:
using System;
using Skyiv.Extensions;
class Tester
{
static void Main()
{
try
{
foreach (var x in new decimal[] {
-100, -66, -65, -1, 0, 1, 2.5m, 16, 66.5421m, 67 })
Console.WriteLine("{0,-30}: exp({1})", x.Exp(), x);
}
catch (Exception ex) { Console.WriteLine(ex.Message); }
}
}
运行结果如下所示:
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的值。这些值是如何得到的呢?这很简单,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 的自然对数。
分享到:
相关推荐
提出一种基于传感器温度补偿的双指数函数拟合算法,一方面采用双指数函数对非线性的温度系数曲线进行补偿,另一方面在双指数函数拟合的算法中,提出一种具有高精度初值的交替迭代法。该方法首先利用四组数据点计算出...
指数函数和对数函数的高精度快速算法_徐洋 学习交流,禁止商业传播
本文基于FPGA实现三角函数、反三角函数以及指数函数计算,分别采用了cordic算法和切比雪夫逼近算法,比较了迭代次数达到误差精度10^-6. 建立已知角度θ,求解sinθ、cosθ的数学模型。 建立已知弧度θ,求解arctan...
本文在这些算法的基础上提出一种双指数函数模型的温度补偿算法,有以下优点:(1)指数函数具有无限阶的泰勒展开式,因此双指数函数在对诸如传感器温度系数曲线这类非线性曲线的拟合上可以达到很高的精度。...
计算张量指数函数的广义逆张量ε-算法.docx
微带结构空域闭式格林函数算法研究,耿海贤,高正平,离散复镜像法快速计算微带结构的空域闭式格林函数,采用广义函数束方法将谱域格林函数表示成指数函数的线性组合,并使用Sommerfeld��
《特殊函数计算手册》较系统地阐述了各种特殊函数的定义、数学性质、算法、数表和程序。由特定微分方程的解定义的特殊函数有正交多项式(如Chebyshev、Laguerre和Hermite多项式),Gamma函数,Legendre函数类,...
【摘要】文章首先介绍CORDIC算法双曲系统的基本原理及其计算模式,对COEDIC内核及前处理单元做了详细分析。 【关键词】坐标旋转数字诗算方法;指数函数;流水线
适合计算李雅普诺夫指数,经典沃夫算法的李雅普诺夫指数
利用遗传算法优化oif指数进行高光谱波段选择 最佳指数因子 高光谱波段选择
正弦余弦算法是一种新型智能优化算法,利用正弦...结果表明指数函数递减的正弦余弦算法具有更高的计算精度和更快的收敛速度。最后以协同过滤推荐算法中相似度函数的计算为应用对象,进一步验证新算法的可行性和有效性。
构建了引导进化算法收敛的指数函数曲线模型,给出了模型的参数计算方法.设计了一种具有全局变异和局部变异算子的进化模糊聚类算法,根据全局变异前后个体适应度值和分量值的变化趋势,实现定向变异,并给出了算法的种群...
文中提出了基于基尼指数的分布估计算法,采用实数编码直接对连续随机变量建模,并引入了基尼指数,设计了可以随着进化代数的变化动态调整子代种群的扰动因子函数.实验结果表明,该算法与其他同类算法相比优化精度有了...
带入数据可以直接运行 可以求出log-log的双对数函数关系图
的时间是关于问题规模的指数函数存在指数爆炸的问题。解决 TSP 问 题我们最直观的想法就是遍历整个图找出所有的 Hamilton 回路再进行 比较、寻优。对于一个具有 n 个顶点的对称完全图而言要从 2)!1(− n 个 ...
基于Putzer算法,将n维空间广义矩阵指数函数的计算转化为一维空间对应的齐次动力学方程求解问题.结果表明,该方法降低了空间维数及计算难度,得到了广义矩阵指数函数的显示表达及算法.
该函数用来计算时间序列的最大Lyapunov 指数--Wolf 方法
资源中包含3个函数(hurst,lors,RSana),可分别计算该指数。 HURST: Calculate the Hurst exponent using R/S analysis lors:Function to compute Lo's (1991) modified rescaled range statistic RSana:...
matlab复变函数指数函数代码set-mifs 基于集合互信息的Matlab特征选择算法的Matlab实现 介绍 在文献中已经提出了使用互信息(MI)来确定模式识别任务中特征的显着性的思想的许多变体。 但是,它们有其局限性:在变量...
CART算法的全称是分类回归树算法,他是一个二元分类,采用的是类似于熵的基尼指数作为分类决策,形成决策树后之后还要进行剪枝,我自己在实现整个算法的时候采用的是代价复杂度算法,详细介绍链接 KNN K最近邻算法...