这个问题是21icbbs 上的一个网友提出的，我第一反应就是迭代。在很多场合下，比如计算热电流thermal current，或是其他一些使用如下公式的应用：
y(t) = K * (1 - exp(-t / T))) (1)
y(t) = K * exp(-t / T)) (2)
where T is the time constant of the system and K is the magnitude.
公式（1），(2)不能直接应用在实际计算中，因为 t并没有确切的数值。 对 (1) 微分后得：
dy = - K * exp(-t / T)) /T = ( K - K*(1 – exp(-t/T))) / T
= (K – y) / T (3)
对 (2) 微分后得：
dy = - K * exp(-t / T)) /T = - y /T (4)
当 K=0, 则 （3）成为 （4）。假设采样时间为 Ts, 把（3） 转成迭代方程:
y(k+1) = y(k) + Ts/T * (K – y(k)) (5)
只要知道初始值 y(0) 就可计算，与时间 t 无关。如果在 s domain ，转换则更简单：把 s = (1- z)/Ts 或 s = 2/Ts * (1+z)/(1-z) 带入下列公式，可得迭代公式：
y(s) = K/(s + 1/T) * x(s) (6)
再看看 exp(x) 公式，如果x 是16-bit整数（网友的要求），那么可以构造出一个矩阵，把 exp(x) 变成乘法问题。
float Coeff[4][16] = {
{ 1, exp(0x0001), exp(0x0002), …., exp(0x000F) },
{ 1, exp(0x0010), exp(0x0020), …., exp(0x00F0) },
{ 1, exp(0x0100), exp(0x0200), …., exp(0x0F00) },
{ 1, exp(0x1000), exp(0x2000), …., exp(0xF000) }
};
float exp(uint16 x)
{
return Coeff[0][x & 0x0F] *
Coeff[1][(x>>4) & 0x0F] *
Coeff[2][(x>>8) & 0x0F] *
Coeff[3][(x>>12) & 0x0F];
}
经过一些网友的提示，我看了一下 CORDIC 算法，上面的方法是把exp(x) 变成乘法，如果能精心挑选x分解项，把 Coeff 变成 2^N，那么意味着在 FPGA 以及 DSP 上只要使用加减和左移右移就能达到目的，这就是CORDIC 的精髓。X 可分解为
x = k0 + k1 + … + kn ,
exp(x) = exp(k0) * exp(k1) * … * exp(kn)
= 2^N0 * 2^N1 * … * 2 * (1+ ½) * (1+ 1/4) * (1+ 1/8) …
注意当xn 趋向0 时，exp(xn) 趋向1, 所以当 0 < x < ln(2) 时，无法单纯的左移右移, 这时可以让 exp(kn) = 1 + 1/(2^N), 则运算变成
y*(1 + 1/(2^N)) = y + (y >> N)
验证程序如下：
#include "stdafx.h"
#include "math.h"
#define LN2 0.69314718055994530941723212145818f
float C1[] =
{
128*LN2, 64*LN2, 32*LN2, 16*LN2, 8*LN2, 4*LN2, 2*LN2, LN2
};
float C2[] =
{
0.405465108f, 0.223143551f, 0.117783036f, 0.060624622f,
0.030771659f, 0.015504187f, 0.00778214f, 0.00389864f,
0.00195122f, 0.000976086f, 0.000488162f, 0.000244111f
};
float Exp(float x)
{
float y = 1.0;
int i;
for (i=0; i<8; i++) {
if (x >= C1[i]) {
x -= C1[i];
y *= 1L <<(1<<(7-i));
}
}
for (i=0; i<12; i++) {
if (x >= C2[i]) {
x -= C2[i];
y += y / (1<<(i+1));
}
}
return y;
}
int _tmain(int argc, _TCHAR* argv[])
{
float y1 = Exp(5.0);
float y2 = (float) exp(5.0);
printf("%4.5f %4.5f", y1, y2);
return 0;
}
用定点计算的程序
#include "stdafx.h"
#include "math.h"
#define uint32 unsigned int
#define int32 int
#define ONE 65536
uint32 C1[] =
{
726817, 363408, 181704, 90852, 45426
};
uint32 C2[] =
{
26572, 14623, 7719, 3973, 2016, 1016,
510, 256, 128, 64, 32, 16, 8, 4, 2
};
// Convert: uint32 = float * 65536
// ONE = 0x00010000 <---> 1.0 in float
// ln(65536) = 11.090,
// so x must be less than 11.090*ONE = 726817
uint32 Exp(uint32 x)
{
uint32 y = ONE;
if (x >= C1[0]) { x -= C1[0]; y <<= 16; }
if (x >= C1[1]) { x -= C1[1]; y <<= 8; }
if (x >= C1[2]) { x -= C1[2]; y <<= 4; }
if (x >= C1[3]) { x -= C1[3]; y <<= 2; }
if (x >= C1[4]) { x -= C1[4]; y <<= 1; }
for (int i=0; i<16; i++) {
if (x >= C2[i]) { x -= C2[i]; y += (y >> (i+1)); }
}
return y;
}
int _tmain(int argc, _TCHAR* argv[])
{
float y1 = (float) Exp((uint32) 5.0f * ONE) / ONE;
float y2 = (float) exp(5.0);
printf("%4.5f %4.5f", y1, y2);
return 0;
}
结果：
148.40210 148.41316
148.41339 148.41316
{var%20f='http://v.t.sina.com.cn/share/share.php?appkey=1515056452',u=z||d.location,p=['&url=',e(u),'&title=',e(t||d.title),'&source=',e(r),'&sourceUrl=',e(l),'&content=',c||'gb2312','&pic=',e(p||'')].join('');function%20a(){if(!window.open([f,p].join(''),'mb',['toolbar=0,status=0,resizable=1,width=440,height=430,left=',(s.width-440)/2,',top=',(s.height-430)/2].join('')))u.href=[f,p].join('');};if(/Firefox/.test(navigator.userAgent))setTimeout(a,0);else%20a();})(screen,document,encodeURIComponent,'','','https://www.xiaopingtou.cn//data/attach/topic/topicKPo7gB.jpg', '推荐 嵌入式小菜 的文章《exp 函数的数值计算方法》','https://www.xiaopingtou.net/article-80510.html','页面编码gb2312|utf-8默认gb2312'));){kind=link}