6

我正在尝试确定我的一种算法的渐近运行时间，该算法使用指数，但我不确定如何以编程方式计算指数。

我专门寻找用于双精度浮点数的 pow() 算法。

4

6 回答 6

14

我有机会查看 fdlibm 的实现。评论描述了使用的算法：

 *                    n
 * Method:  Let x =  2   * (1+f)
 *      1. Compute and return log2(x) in two pieces:
 *              log2(x) = w1 + w2,
 *         where w1 has 53-24 = 29 bit trailing zeros.
 *      2. Perform y*log2(x) = n+y' by simulating muti-precision
 *         arithmetic, where |y'|<=0.5.
 *      3. Return x**y = 2**n*exp(y'*log2)

然后是处理的所有特殊情况的列表（0、1、inf、nan）。

在所有特殊情况处理之后，代码中最密集的部分涉及log2和2**计算。并且其中任何一个都没有循环。因此，尽管浮点原语很复杂，但它看起来像是一个渐近恒定时间算法。

欢迎浮点专家（我不是其中之一）发表评论。:-)

于 2008-10-03T00:16:34.000 回答
3

除非他们发现了更好的方法，否则我相信三角函数、对数函数和指数函数（例如指数增长和衰减）的近似值通常使用算术规则和泰勒级数展开计算，以产生准确的近似结果在要求的精度范围内。（有关幂级数、泰勒级数和麦克劳林级数展开函数的详细信息，请参阅任何微积分书籍。）请注意，我已经有一段时间没有这样做了，所以我不能告诉你，例如，确切地如何计算您需要包含的系列中的项数可以保证误差小到在双精度计算中可以忽略不计。

例如，e^x 的 Taylor/Maclaurin 级数展开是这样的：

      +inf [ x^k ]           x^2    x^3      x^4        x^5
e^x = SUM  [ --- ] = 1 + x + --- + ----- + ------- + --------- + ....
      k=0  [  k! ]           2*1   3*2*1   4*3*2*1   5*4*3*2*1

如果您采用所有项（k 从 0 到无穷大），则此展开是准确且完整的（无错误）。

但是，如果您没有将所有项都取到无穷大，而是在说 5 个项或 50 个项或其他任何项之后停止，您会产生一个与实际 e^x 函数值相差一个余数的近似结果，这很容易计算。

指数的好消息是它很好地收敛并且它的多项式展开的项很容易迭代编码，所以你可能（重复，可能- 记住，已经有一段时间了）甚至不需要预先计算你有多少项需要保证您的错误小于精度，因为您可以在每次迭代时测试贡献的大小，并在它变得足够接近零时停止。在实践中，我不知道这个策略是否可行——我必须尝试一下。有些重要的细节我早就忘记了。诸如：机器精度，机器误差和舍入误差等。

另外，请注意，如果您没有使用 e^x，但您正在使用另一个基数（如 2^x 或 10^x）进行增长/衰减，则近似多项式函数会发生变化。

于 2008-10-03T01:33:32.583 回答
1

对于整数指数，通常的方法是将 a 提高到 b，如下所示：

result = 1
while b > 0
  if b is odd
    result *= a
    b -= 1
  b /= 2
  a = a * a

它的指数大小通常是对数的。该算法基于不变量“a^b*result = a0^b0”，其中 a0 和 b0 是 a 和 b 的初始值。

对于负或非整数指数，需要对数和近似值以及数值分析。运行时间将取决于使用的算法和库的调整精度。

编辑：由于似乎有一些兴趣，这里有一个没有额外乘法的版本。

result = 1
while b > 0
  while b is even
    a = a * a
    b = b / 2
  result = result * a
  b = b - 1
于 2008-10-02T23:01:06.443 回答
1

您可以使用 exp(n*ln(x)) 来计算 x n。x 和 n 都可以是双精度浮点数。自然对数和指数函数可以使用泰勒级数计算。在这里你可以找到公式：http ://en.wikipedia.org/wiki/Taylor_series

于 2008-10-07T11:42:09.703 回答
0

如果我正在编写一个针对 Intel 的 pow 函数，我将返回 exp2(log2(x) * y)。英特尔的 log2 微代码肯定比我能编写的任何代码都快，即使我还记得我第一年的微积分和研究生院的数值分析。

于 2008-10-03T00:30:26.900 回答
0

e^x = (1 + 分数) * (2^exponent), 1 <= 1 + 分数 < 2

x * log2(e) = log2(1 + 分数) + 指数，0 <= log2(1 + 分数) < 1

指数 = floor(x * log2(e))

1 + 分数 = 2^(x * log2(e) - 指数) = e^((x * log2(e) - 指数) * ln2) = e^(x - 指数 * ln2), 0 <= x - 指数* ln2 < ln2

于 2021-08-27T21:36:21.907 回答