近似 e^x

Posted 2023-02-18

【中文标题】近似 e^x

【英文标题】：Approximate e^x

【发布时间】：2011-10-22 12:04:12

【问题描述】：

我想近似 e^x 函数。

是否可以使用基于多个样条类型的方法来做到这一点？即在 x₁ 和 x₂ 之间，然后是

y₁ = a₁x + b₁，在x₂之间> 和 x₃,

然后

y₂ = a₂x + b₂

等

这是用于专用 fpga 硬件而不是通用 CPU。因此，我需要自己创建函数。准确度就不是那么重要了。此外，我真的买不起超过一个乘法电路和/或多个移位/加法器。我也想要比 CORDIC 函数小得多的东西，实际上大小很关键。

【问题讨论】：

您打算在多大范围内逼近 x 值？

默认答案：power series

您在 C++ 标准中有 exp() 函数。你为什么避免使用它？通常它有很好的速度。

递归近似不适合我的应用。潜在的最大范围是 0-4095，但可以缩放到更小的值。我的直觉是我需要大约 4 到 6 位的精度

我的应用程序实际上不是 C 或 C++，它的专用硬件，所以我自己滚动函数。电源功能很好，但我更喜欢操作较少的东西。

【参考方案1】：

像这样使用公式的策略怎么样

e^x = 2 ^x/ln(2)

预先计算1/ln(2) 将此常数乘以您的论点（1 次乘法） 使用二进制移位将 2 提高到幂的整数部分（假设为 exp+mantissa 格式） 根据 2 的小数次方余数进行调整（可能是二次乘法）

我意识到这不是一个完整的解决方案，但它只需要一次乘法，并将剩余的问题减少到近似 2 的分数幂，这应该更容易在硬件中实现。

此外，如果您的应用程序足够专业，您可以尝试重新导出将在您的硬件上运行的所有数字代码，使其处于一个基本的e 数字系统中并实现您的浮动点硬件也可以在基础 e 中工作。那么根本不需要转换。

【讨论】：

感谢 Lucas - 这非常适合我的需求，甚至比我所希望的还要好。非常感谢！

很高兴听到。听起来你有一些有趣的设计权衡。

@trican 有一篇很好的论文介绍了如何使用查找表和定点算法实现这种身份和范围缩减，以实现单精度浮点的合理精度：loria.fr/~detreyje/publications/DetDin_fpt_2005.pdf

PDF 的替代链接：perso.citi-lab.fr/fdedinec/recherche/publis/2005-FPT.pdf

【参考方案2】：

如果x 是一个整数，您可以一遍又一遍地乘以e。

如果x不是整数，可以使用上述方法计算出e^floor(x)，然后乘以一个小的修正项。可以使用多种近似方法轻松计算该校正项。一种这样的方法是：

e^f ≈ 1 + f(1 + f/2(1 + f/3(1 + f/4)))，其中 f 是 x 的小数部分

这来自 e^x 的（优化的）幂级数展开式，这对于 x 的小值非常准确。如果您需要更高的准确性，只需在该系列中添加更多术语。

这个math.stackexchange 问题包含一些额外的聪明答案。

编辑：请注意，计算 eⁿ 有一种更快的方法，称为 exponentiation by squaring。

【讨论】：

整数解的最佳解不是这个 O(n) 解。分而治之的算法（预）计算 e^1、e^2、e^4、e^8 等。然后您取与x 中的位相对应的因子。这是 O(logN)。 IE。对于 x=255，这只需要 8 次乘法而不是 254 次。

谢谢 - 但我希望尽量减少乘法运算，我只想要一个乘法运算

但是为什么？您是否实际上看到了性能问题，或者这是过早的优化？

@Jonathan - 它不是用于 CPU，而是用于专用硬件。我已经更新了上面的问题以澄清这一点。很抱歉造成混乱

@Jonathan 因为拥有 O(n) 指数函数显然会导致性能不佳。过早的优化在系统层面上还不错。

【参考方案3】：

首先，是什么激发了这种近似？换句话说，直截了当的exp(x)到底有什么问题？

也就是说，exp(x) 的典型实现是

找到一个整数k 和浮点数r，使得x=k*log(2) + r 和r 在-0.5*log(2) 和0.5*log(2) 之间。 通过这种减少，exp(x) 为 2^k*exp(r)。 计算 2^k 很简单。 exp(x) 的标准实现使用 Remes 类型的算法得出一个近似于 exp(r) 的极小极大多项式。 您也可以这样做，但使用降阶多项式。

这里是关键：无论你做什么，你的函数都会比调用exp() 慢很多，慢得多。 exp() 的大部分功能都在您计算机的数学协处理器中实现。在软件中重新实现该功能，即使精度降低，也将比仅使用 exp() 慢一个数量级。

【讨论】：

Remez* 和大多数实际使用以边界为中心的 Pade 近似值，以便在此范围内的误差尽可能小。给定输入x 的误差等于有界误差乘以2^k，当输入很大时，这通常会破坏这些近似值...输入减去反函数的迭代改进求根方法。

为什么r 应该位于-0.5log(2) 和0.5log(2) 而不是(0, 1) 之间？

【参考方案4】：

对于硬件，如果您需要精确到位，我可以为您提供一个很棒的解决方案。 （否则只需像上面那样做一个近似值）。恒等式是 exp(x) = cosh(x) + sinh(x)，即双曲正弦和余弦。问题是双曲正弦和余弦可以使用 CORIC 技术计算，最重要的是，它们是 FAST CORDIC 函数之一，这意味着它们看起来几乎像乘法而不是像除法！

这意味着对于数组乘法器的面积，您可以在 2 个周期内计算任意精度的指数！

查找 CORDIC 方法 - 这对于硬件实现来说非常棒。

另一种硬件方法是使用小表格和其他人提到的公式：exp(x + y) = exp(x) * exp(y)。您可以将数字分解为小的位域（例如一次 4 位或 8 位），然后只需查找该位域的指数即可。可能只对狭窄的计算有效，但这是另一种方法。

【参考方案5】：

http://martin.ankerl.com/2007/02/11/optimized-exponential-functions-for-java/ 使用 Schraudolph 的方法 (http://nic.schraudolph.org/pubs/Schraudolph99.pdf) 在 Java 中：

public static double exp(double val) 
    final long tmp = (long) (1512775 * val) + (1072693248 - 60801);
    return Double.longBitsToDouble(tmp << 32);

和 https://math.stackexchange.com/a/56064（寻找 Pade 近似值）。

【讨论】：

感谢@jdberton 添加此内容和链接。这个方法看起来很有趣，但是你确定上面的代码 sn-p 是正确的吗？我尝试了一些值，结果似乎并不接近？

我认为它对于大值是不准确的。您可能会通过一些工作找到更好的 Pade 近似值以获得更好的范围。它对我有用，因为我不需要任何确切的东西。

Schraudolphs 方法非常完美。如果准确性可以接受，我认为它不会变得更快。在他的论文中，他确定平均相对误差约为 4%。来源：nic.schraudolph.org/pubs/Schraudolph99.pdf

这里是 Schraudolph 方法的更现代的实现，使用单点浮点而不是双精度（这是一种浪费，因为只写入双精度的高 32 位）。 machinedlearnings.com/2011/06/…

【参考方案6】：

这不是您要求的平滑样条插值，而是计算效率：

float expf_fast(float x) 
   union  float f; int i;  y;
   y.i = (int)(x * 0xB5645F + 0x3F7893F5);
   return (y.f);

绘图输出

【参考方案7】：

Wolfram 提出了一些根据级数等近似它的好方法：

Wolfram page for e^x

Taylor Series 上的***页面还显示了 e^x 在 0 附近扩展的示例：

【讨论】：

"替代表示：e^x=z^x for e=z" :D

【参考方案8】：

或者您可以在 C 中执行 pow(M_E, x)。（某些平台没有定义 M_E；在这些平台上，您可能必须手动指定 e 的值，大约是2.71828182845904523536028747135266249775724709369995.)

（正如 David 在 cmets 中指出的那样，exp(x) 将比 pow(M_E, x) 更有效。同样，大脑还没有打开。）

您是否有一个用例，其中 e^x 的计算是一个已证实的瓶颈？如果没有，您应该首先编写可读性代码；只有在明显的方法太慢时才尝试这些优化。

【讨论】：

pow(M_E, x)?严重地？ pow(a,b) 通常实现为 exp(b*log(a))。使用pow 是一个减速带，而不是加速。

这就是我的观点——先正确编写代码，然后看看它的性能。原始问题中没有任何地方说这被称为每秒一百万次或类似的东西，因此性能是否会成为问题并不是很明显。

无论性能如何，exp(x) 都是比pow(M_E, x) 更简单（且更便携！）的解决方案。即使pow() 更快，使用它而不是exp() 也将是过早的优化。

非常正确，我已经更新了我的答案以反映大卫的更正。你能说我还没有喝足够的咖啡吗？ :)

【参考方案9】：

当然这是“可能的”。有几个问题。

您对精度有什么要求？

您愿意使用高阶样条吗？

您愿意为此花费多少内存？足够小区间内的线性函数可以将指数函数逼近到所需的任何精度，但它可能需要非常小的区间。

编辑：

鉴于提供的其他信息，我进行了快速测试。范围缩小总是可以用于指数函数。因此，如果我希望计算任意 x 的 exp(x)，那么我可以将问题改写为...

y = exp(xi + xf) = exp(xi)*exp(xf)

其中 xi 是 x 的整数部分，xf 是小数部分。整数部分很简单。以二进制形式计算 xi，然后重复平方和乘法允许您在相对较少的操作中计算 exp(xi)。 （其他技巧，使用 2 的幂和其他间隔可以为您提供更快的速度。）

现在剩下的就是计算 exp(xf)。我们可以使用带有线性段的样条来计算 exp(xf)，在区间 [0,1] 上只有 4 个线性段，精度为 0.005？

最后一个问题由我几年前编写的一个函数解决，该函数将使用给定阶数的样条曲线逼近一个函数，使其在最大误差的固定容差范围内。此代码需要区间 [0,1] 上的 8 段，以通过分段线性样条函数实现所需的容差。如果我选择将区间进一步减小到 [0,0.5]，我现在可以达到规定的容差。

所以答案很简单。如果您愿意进行范围缩小以将 x 减小到区间 [0.0.5]，然后进行适当的计算，那么您可以通过 4 段的线性样条曲线达到要求的精度。

最后，尽管使用硬编码的指数函数总是会更好。上面提到的所有操作肯定会比你的编译器提供的要慢，如果 exp(x) 可用。

【讨论】：

非常感谢您的详细回复。进一步思考，我可以容忍更高的误差范围，可能高达 0.05，甚至可能是 0.1。我之前曾将样条曲线与范围缩减一起用于其他功能，但在这种情况下，我认为 Lucas 的上述回答更适合较低的精度要求。此外，关键点是在硬件“编译器”中没有直接实现指数函数。即我没有在 CPU 上工作

【参考方案10】：

这不适用于定制 FPGA，但值得一提。

http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html

以及源代码：

https://code.google.com/archive/p/fastapprox/downloads

“更快”的实现只涉及 3 个步骤（乘法、加法、将 float 转换为 int）和最终转换回 float。根据我的经验，它的准确率为 2%，如果您不关心实际值但在对数似然最大化迭代中使用该值，这可能就足够了。