pow计算

ab=elnab=eb∗lnaa^b=e^{ln\ a^b}=e^{b*ln\ a} ab=eln ab=eb∗ln a所以任意指数都可以有自然指数和对数来求得。

exp的算法

形如y=exy=e^xy=ex的指数函数用底层实现，最直接的是泰勒逼近，然后找到一篇国内的paper指数函数e^x的快速计算方法，分析下来感觉不错，记录一下

泰勒级数展开

exe^xex的泰勒展开为：ex=1+x1!+x22!+x33!+x44！+....e^x=1+\frac{x}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4！}+.... ex=1+1!x+2!x2+3!x3+4！x4+....分析下来，如果∣x∣<1|x|<1∣x∣<1，高次方很快收敛，那么相对容易计算这个逼近值。但一般情况比较难确定。

快速算法

y=exy=e^xy=ex这个等式两端取2为底的对数log2y=xlog2e=ax,a=log2elog_2\ y=xlog_2\ e=ax, a=log_2\ e log2 y=xlog2 e=ax,a=log2 e则经过变换y=2axy=2^{ax}y=2ax怎么理解这个公式呢？即将自然数底的指数变成呢二进制底的指数，因为a=log2ea=log_2\ ea=log2 e是常数，好像结合计算机的二进制方法可以算出来呢。分析假设axaxax的整数部分为nnn，小数部分为E=ax−nE=ax-nE=ax−n，这样重写y=2n+E=2n.2Ey=2^{n+E}=2^n.2^Ey=2n+E=2n.2E，那么只要计算出axaxax保留整数部分，剩下的小数部分可以变换回自然指数展开的形式，求得后与整数部分相乘即可。小数部分我们设为p=2E=eln2E=eEln2p=2^E=e^{ln2^E}=e^{Eln2}p=2E=eln2E=eEln2，假设x00=Eln2x_{00}=Eln2x00=Eln2因为E最大值小于1，所以0<x00<ln2=0.69310<x_{00}<ln2=0.69310<x00<ln2=0.6931这部分可以利用泰勒展开来逼近呢，但是由于结果和整数部分相乘，所以误差计算的时候不能忽略。

进一步减少尾数

把x00x_{00}x00进一步分解，可得如下表达：
p=eEln2=ex00=ex0+δxp=e^{Eln2}=e^{x_{00}}=e^{x_0+\delta x} p=eEln2=ex00=ex0+δx遵循原文中的表达，采用3位16进制小数进行分解，则x00=[0.p1p2p3+0.000δ]hexx_{00}=[0.p_1p_2p_3+0.000\delta]_{hex} x00=[0.p1p2p3+0.000δ]hexpxp_xpx可以将x_{00}以此左移4位保留整数部分求得，其中p0<ln2∗16=0.6931∗16p_0<ln2 \ *16=0.6931*16p0<ln2 ∗16=0.6931∗16。
剩下的eδxe^{\delta x}eδx就可以用低阶的泰勒展开求解，然后将整数部分，e0.pxe^{0.p_x}e0.px部分和eδxe^{\delta x}eδx统统乘起来就是所得的结果。

利用双精度表达提高运算效率

双精度需要8bytes保存数位，最高位符号位S，接下来11bit存储2底幂指数n，生下来的52位Mask存储2底小数（<2)的部分。转换公式为：decimal=(−1)S∗2n−1023∗(1+∑i=152Mi∗2−i)=(−1)S∗2n−1023∗(∑i=052Mi∗2−i)decimal=(-1)^S*2^{n-1023}*(1+\sum_{i=1}^{52}M_i*2^{-i})\\=(-1)^S*2^{n-1023}*(\sum_{i=0}^{52}M_i*2^{-i}) decimal=(−1)S∗2n−1023∗(1+i=1∑52Mi∗2−i)=(−1)S∗2n−1023∗(i=0∑52Mi∗2−i)
咋看起来和上面公式有些出入，但实际上y=2n+E=2n.2E2E=(∑i=052Mi∗2−i)y=2^{n+E}=2^n.2^E\\2^E=(\sum_{i=0}^{52}M_i*2^{-i})y=2n+E=2n.2E2E=(i=0∑52Mi∗2−i)完美融合，因为202^020是不保存在双精度浮点数位里面的。因为幂指数没有负数，S位填0，即：

nnn位减去1023填入幂指数11bit；
将e0.pxe^{0.p_x}e0.px部分和eδxe^{\delta x}eδx求解出来得积换算成2E2^E2E的形式（∗252*2^{52}∗252）填入小数部分。

所以算好的数据填入一个双精度的联合体，答案就出来了。

实验

利用自己设计的3*16表格，作为输入，简单测试了精确度，和C的库保持的很好。

x=1.00000000, expw=2.7182818285,  exp(x)=2.7182818285
x=1.06449446, expw=2.8993728614,  exp(x)=2.8993728614
x=1.13314845, expw=3.1054183887,  exp(x)=3.1054183887
x=1.20623025, expw=3.3408666501,  exp(x)=3.3408666501
x=1.28402542, expw=3.6111468782,  exp(x)=3.6111468782
x=1.36683794, expw=3.9229265384,  exp(x)=3.9229265384
x=1.45499141, expw=4.2844466818,  exp(x)=4.2844466818
x=1.54883030, expw=4.7059623913,  exp(x)=4.7059623913
x=1.64872127, expw=5.2003257648,  exp(x)=5.2003257648
x=1.75505466, expw=5.7837638563,  exp(x)=5.7837638563
x=1.86824596, expw=6.4769256265,  exp(x)=6.4769256265
x=1.98873747, expw=7.3063035064,  exp(x)=7.3063035064
x=2.11700002, expw=8.3061816659,  exp(x)=8.3061816659
x=2.25353479, expw=9.5213323069,  exp(x)=9.5213323069
x=2.39887529, expw=11.0107855170,  exp(x)=11.0107855170
x=2.55358946, expw=12.8531569481,  exp(x)=12.8531569481
x=1.00000000, expw=2.7182818285,  exp(x)=2.7182818285
x=1.00391389, expw=2.7289417300,  exp(x)=2.7289417300
x=1.00784310, expw=2.7396854025,  exp(x)=2.7396854025
x=1.01178768, expw=2.7505136706,  exp(x)=2.7505136706
x=1.01574771, expw=2.7614273686,  exp(x)=2.7614273686
x=1.01972323, expw=2.7724273406,  exp(x)=2.7724273406
x=1.02371432, expw=2.7835144407,  exp(x)=2.7835144407
x=1.02772102, expw=2.7946895334,  exp(x)=2.7946895334
x=1.03174341, expw=2.8059534932,  exp(x)=2.8059534932
x=1.03578154, expw=2.8173072053,  exp(x)=2.8173072053
x=1.03983547, expw=2.8287515653,  exp(x)=2.8287515653
x=1.04390527, expw=2.8402874797,  exp(x)=2.8402874797
x=1.04799100, expw=2.8519158657,  exp(x)=2.8519158657
x=1.05209272, expw=2.8636376517,  exp(x)=2.8636376517
x=1.05621050, expw=2.8754537771,  exp(x)=2.8754537771
x=1.06034439, expw=2.8873651929,  exp(x)=2.8873651929
x=1.00000000, expw=2.7182818285,  exp(x)=2.7182818285
x=1.00024417, expw=2.7189456335,  exp(x)=2.7189456335
x=1.00048840, expw=2.7196097629,  exp(x)=2.7196097629
x=1.00073269, expw=2.7202742166,  exp(x)=2.7202742166
x=1.00097704, expw=2.7209389950,  exp(x)=2.7209389950
x=1.00122145, expw=2.7216040983,  exp(x)=2.7216040983
x=1.00146592, expw=2.7222695265,  exp(x)=2.7222695265
x=1.00171045, expw=2.7229352800,  exp(x)=2.7229352800
x=1.00195503, expw=2.7236013590,  exp(x)=2.7236013590
x=1.00219968, expw=2.7242677635,  exp(x)=2.7242677635
x=1.00244439, expw=2.7249344940,  exp(x)=2.7249344940
x=1.00268916, expw=2.7256015504,  exp(x)=2.7256015504
x=1.00293398, expw=2.7262689330,  exp(x)=2.7262689330
x=1.00317887, expw=2.7269366421,  exp(x)=2.7269366421
x=1.00342382, expw=2.7276046778,  exp(x)=2.7276046778
x=1.00366882, expw=2.7282730404,  exp(x)=2.7282730404

ln(x)计算

作为exp(x)exp(x)exp(x)的逆运算，上面的步骤可以参考，先看一下关于lnlnln的泰勒展开：
ln(1+x)=∑0∞(−1)nn+1xn+1ln(1+x)=\sum_0^\infty\frac{(-1)^n}{n+1}x^{n+1} ln(1+x)=0∑∞n+1(−1)nxn+1回顾decimal=(−1)S∗2n−1023∗(1+∑i=152Mi∗2−i)decimal=(-1)^S*2^{n-1023}*(1+\sum_{i=1}^{52}M_i*2^{-i}) decimal=(−1)S∗2n−1023∗(1+i=1∑52Mi∗2−i)，抛开符号位（对数的输入都是正的）对这样一个对数运算，ln(decimal)=ln(2n−1023∗(1+∑i=152Mi∗2−i))=ln2n−1023+ln(1+∑i=152Mi∗2−i)=(n−1023)∗ln2+ln(1+∑i52Mi∗2−i)ln(decimal)=ln(2^{n-1023}*(1+\sum_{i=1}^{52}M_i*2^{-i}))\\=ln\ 2^{n-1023}+ln(1+\sum_{i=1}^{52}M_i*2^{-i})\\=(n-1023)*ln\ 2+ln(1+\sum_i^{52}M_i*2^{-i}) ln(decimal)=ln(2n−1023∗(1+i=1∑52Mi∗2−i))=ln 2n−1023+ln(1+i=1∑52Mi∗2−i)=(n−1023)∗ln 2+ln(1+i∑52Mi∗2−i)上式的第一项是输入乘以常数项ln2ln\ 2ln 2，ln(1+∑i=152Mi∗2−i)ln(1+\sum_{i=1}^{52}M_i*2^{-i})ln(1+∑i=152Mi∗2−i)作为一个整体用泰勒级数求解，结果可得。

实验

因为泰勒展开是正负交替，所以收敛性能不如单边稳定，以100级泰勒展开为例，对比math.h库的结果如下：

x=0.10000000, loge=-2.3025850930,  log(x)=-2.3025850930
x=0.20000000, loge=-1.6094379124,  log(x)=-1.6094379124
x=0.30000000, loge=-1.2039728043,  log(x)=-1.2039728043
x=0.40000000, loge=-0.9162907319,  log(x)=-0.9162907319
x=0.50000000, loge=-0.6931471806,  log(x)=-0.6931471806
x=0.60000000, loge=-0.5108256238,  log(x)=-0.5108256238
x=0.70000000, loge=-0.3566749439,  log(x)=-0.3566749439
x=0.80000000, loge=-0.2231435513,  log(x)=-0.2231435513
x=0.90000000, loge=-0.1053605157,  log(x)=-0.1053605157
x=1.00000000, loge=-0.0049750012,  log(x)=-0.0000000000
x=1.10000000, loge=0.0953101798,  log(x)=0.0953101798
x=1.20000000, loge=0.1823215568,  log(x)=0.1823215568
x=1.30000000, loge=0.2623642645,  log(x)=0.2623642645
x=1.40000000, loge=0.3364722366,  log(x)=0.3364722366
x=1.50000000, loge=0.4054651081,  log(x)=0.4054651081
x=1.60000000, loge=0.4700036292,  log(x)=0.4700036292
x=1.70000000, loge=0.5306282511,  log(x)=0.5306282511
x=1.80000000, loge=0.5877866649,  log(x)=0.5877866649
x=1.90000000, loge=0.6418537610,  log(x)=0.6418538862
x=2.00000000, loge=0.6931471806,  log(x)=0.6931471806
x=2.10000000, loge=0.7419373447,  log(x)=0.7419373447
x=2.20000000, loge=0.7884573604,  log(x)=0.7884573604
x=2.30000000, loge=0.8329091229,  log(x)=0.8329091229
x=2.40000000, loge=0.8754687374,  log(x)=0.8754687374
x=2.50000000, loge=0.9162907319,  log(x)=0.9162907319
x=2.60000000, loge=0.9555114450,  log(x)=0.9555114450
x=2.70000000, loge=0.9932517730,  log(x)=0.9932517730
x=2.80000000, loge=1.0296194172,  log(x)=1.0296194172
x=2.90000000, loge=1.0647107370,  log(x)=1.0647107370
x=3.00000000, loge=1.0986122887,  log(x)=1.0986122887
x=3.10000000, loge=1.1314021115,  log(x)=1.1314021115
x=3.20000000, loge=1.1631508098,  log(x)=1.1631508098
x=3.30000000, loge=1.1939224685,  log(x)=1.1939224685
x=3.40000000, loge=1.2237754316,  log(x)=1.2237754316
x=3.50000000, loge=1.2527629685,  log(x)=1.2527629685
x=3.60000000, loge=1.2809338455,  log(x)=1.2809338455
x=3.70000000, loge=1.3083328193,  log(x)=1.3083328197
x=3.80000000, loge=1.3350009416,  log(x)=1.3350010667
x=3.90000000, loge=1.3609478574,  log(x)=1.3609765531
x=4.00000000, loge=1.3862943611,  log(x)=1.3862943611
x=4.10000000, loge=1.4109869737,  log(x)=1.4109869737
x=4.20000000, loge=1.4350845253,  log(x)=1.4350845253
x=4.30000000, loge=1.4586150227,  log(x)=1.4586150227
x=4.40000000, loge=1.4816045409,  log(x)=1.4816045409
x=4.50000000, loge=1.5040773968,  log(x)=1.5040773968
x=4.60000000, loge=1.5260563035,  log(x)=1.5260563035
x=4.70000000, loge=1.5475625087,  log(x)=1.5475625087
x=4.80000000, loge=1.5686159179,  log(x)=1.5686159179
x=4.90000000, loge=1.5892352051,  log(x)=1.5892352051
x=5.00000000, loge=1.6094379124,  log(x)=1.6094379124

参考

Numerical Approximations
how is log(x) calculated
指数函数e^x的快速计算方法
DSP Trick: Quick-and-Dirty Logarithms
Fastandcorrectly roundedlogarithmsin double-precision
Exponential and Logarithmic Integrals

研究一下exp, ln, pow的数值计算相关推荐

Eclipse - undefined reference to sin - cos - exp - sqrt - pow
Eclipse - undefined reference to sin - cos - exp - sqrt - pow undefined reference to `pow' undefined ...
log函数 oracle power_Excel之数学函数SQRT/MOD/EXP/LN/RAND
本部分主要包括ABS函数.SQRT函数.SIGN函数.MOD函数.POWER.EXP函数.LN函数.LOG函数.LOG10函数.RAND函数.RANDBETWEEN函数.PI函数.SIN函数.COS函 ...
rand()函数100000随机数_Excel之数学函数SQRT/MOD/EXP/LN/RAND
本部分主要包括ABS函数.SQRT函数.SIGN函数.MOD函数.POWER.EXP函数.LN函数.LOG函数.LOG10函数.RAND函数.RANDBETWEEN函数.PI函数.SIN函数.COS函 ...
《Python程序设计与算法基础教程（第二版）》江红余青松全部章节的课后习题，上机实践，课后答案，案例研究
(还在更新中-) 这篇博客花费了我的大量时间和精力,从创作到维护:若认可本篇博客,希望给一个点赞.收藏并且,遇到了什么问题,请在评论区留言,我会及时回复的这本书对Python的知识点的描述很详细, ...
ln x的matlab表示,ln在matlab中怎么表示
实现thln13算法的matlab程序_数学_自然科学_专业资料.clear a... (x) a^x ln x ax logba cos x tan x cot x ... (2) 指数和对数函数指 ...
ln x的matlab表示,matlab中ln怎么表示
ln x ; y 6 ? lg x ; y 7 ? sin x ; y 8 ? arcsin x . 5 [matlab 命令] >> x=2.2; >> y1=x^3;y1= ...
matlab中ln4怎么表示,matlab里ln怎么表示
matlab中图中颜色和形状表示_数学_自然科学_专业资料.颜色字符串有'c', 'm', 'y', 'r', 'g', 'b', 'w',和'k'.分别表示青,红紫,黄,红,绿,白...... 上下 ...
Math.Pow（）是如何在.NET Framework中实现的？
我一直在寻找一种计算b的有效方法(比如a = 2和b = 50 ). 为了开始,我决定看一下Math.Pow()函数的实现. 但在.NET Reflector中 ,我发现的只有: [MethodImp ...
基于材料数值计算大数据的材料辐照机理发现
点击上方蓝字关注我们基于材料数值计算大数据的材料辐照机理发现任帅1,2, 陈丹丹1,2, 储根深1,2, 白鹤1,2, 李慧昭1, 何远杰1, 胡长军1,2 1 北京科技大学计算机与通信工程学院, ...

研究一下exp, ln, pow的数值计算