admin管理员组

文章数量:1794759

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

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

• pow计算
• exp的算法
• • 泰勒级数展开
  • 快速算法
  • 进一步减少尾数
  • 利用双精度表达提高运算效率
  • 实验
• ln(x)计算
• • 实验
• 参考

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.p1​p2​p3​+0.000δ]hex​pxp_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∑52​Mi​∗2−i)=(−1)S∗2n−1023∗(i=0∑52​Mi​∗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∑52​Mi​∗2−i)完美融合，因为202^020是不保存在双精度浮点数位里面的。因为幂指数没有负数，S位填0，即：

1. nnn位减去1023填入幂指数11bit；
2. 将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)n​xn+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∑52​Mi​∗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∑52​Mi​∗2−i))=ln 2n−1023+ln(1+i=1∑52​Mi​∗2−i)=(n−1023)∗ln 2+ln(1+i∑52​Mi​∗2−i)上式的第一项是输入乘以常数项ln2ln\ 2ln 2，ln(1+∑i=152Mi∗2−i)ln(1+\sum_{i=1}^{52}M_i*2^{-i})ln(1+∑i=152​Mi​∗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

本文标签： 研究一下explnpow的数值计算