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

pow计算

$$a^b=e^{ln{a}^b}=e^{b\ast lna}$$所以任意指数都可以有自然指数和对数来求得。

exp的算法

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

泰勒级数展开

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

快速算法

$y=e^x$这个等式两端取2为底的对数$$lo{g}_{2}y=xlo{g}_{2}e=ax,a=lo{g}_{2}e$$则经过变换$$y=2^{ax}$$怎么理解这个公式呢？即将自然数底的指数变成呢二进制底的指数，因为$a=lo{g}_{2}e$是常数，好像结合计算机的二进制方法可以算出来呢。分析假设$ax$的整数部分为$n$，小数部分为$E=ax-n$，这样重写$$y=2^{n+E}=2^n.2^E$$，那么只要计算出$ax$保留整数部分，剩下的小数部分可以变换回自然指数展开的形式，求得后与整数部分相乘即可。小数部分我们设为$$p=2^E=e^{ln{2}^E}=e^{Eln2}$$，假设$$x_{00}=Eln2$$因为E最大值小于1，所以$$0<x_{00}<ln2=0.6931$$这部分可以利用泰勒展开来逼近呢，但是由于结果和整数部分相乘，所以误差计算的时候不能忽略。

进一步减少尾数

把$x_{00}$进一步分解，可得如下表达：
$$p=e^{Eln2}=e^{x_{00}}=e^{x_0+\delta x}$$遵循原文中的表达，采用3位16进制小数进行分解，则$$x_{00}=[0.p_1{p}_2{p}_3+0.000\delta {]}_{hex}$$$p_x$可以将x_{00}以此左移4位保留整数部分求得，其中$p_0<ln2\ast 16=0.6931\ast 16$。
剩下的$e^{\delta x}$就可以用低阶的泰勒展开求解，然后将整数部分，$e^{0.p_x}$部分和$e^{\delta x}$统统乘起来就是所得的结果。

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

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

1. $n$位减去1023填入幂指数11bit；
2. 将$e^{0.p_x}$部分和$e^{\delta x}$求解出来得积换算成$2^E$的形式（$\ast 2^{52}$）填入小数部分。

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

实验

利用自己设计的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)$的逆运算，上面的步骤可以参考，先看一下关于$ln$的泰勒展开：
$$ln(1+x)=\sum _0^{\infty }\frac{(-1{)}^n}{n+1}{x}^{n+1}$$回顾$$decimal=(-1{)}^S\ast 2^{n-1023}\ast (1+\sum _{i=1}^{52}{M}_i\ast 2^{-i})$$，抛开符号位（对数的输入都是正的）对这样一个对数运算，$$\begin{gathered} ln(decimal)=ln(2^{n-1023}\ast (1+\sum _{i=1}^{52}{M}_i\ast 2^{-i})) \\ =ln{2}^{n-1023}+ln(1+\sum _{i=1}^{52}{M}_i\ast 2^{-i}) \\ =(n-1023)\ast ln2+ln(1+\sum _i^{52}{M}_i\ast 2^{-i}) \end{gathered}$$上式的第一项是输入乘以常数项$ln2$，$ln(1+{\sum }_{i=1}^{52}{M}_i\ast 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