当前位置：文江博客话题详情

指数是如何计算的？

发布于 2024-07-06 18:47:59 字数 88 浏览 12 评论 0原文

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

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

原文

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

逆蝶 2024-07-13 18:47:59

我有机会了解 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** 计算。并且其中任何一个都没有循环。因此，尽管浮点原语很复杂，但它看起来像是一个渐进恒定时间算法。

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

I've had a chance to look at fdlibm's implementation. The comments describe the algorithm used:

 *                    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)

followed by a listing of all the special cases handled (0, 1, inf, nan).

The most intense sections of the code, after all the special-case handling, involve the log2 and 2** calculations. And there are no loops in either of those. So, the complexity of floating-point primitives notwithstanding, it looks like a asymptotically constant-time algorithm.

Floating-point experts (of which I'm not one) are welcome to comment. :-)

回复收藏 0 原文

烟花易冷人易散 2024-07-13 18:47:59

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

例如，e^x 的泰勒/麦克劳林级数展开式是这样的：

      +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）（近似多项式）进行增长/衰减功能变化。

Unless they've discovered a better way to do it, I believe that approximate values for trig, logarithmic and exponential functions (for exponential growth and decay, for example) are generally calculated using arithmetic rules and Taylor Series expansions to produce an approximate result accurate to within the requested precision. (See any Calculus book for details on power series, Taylor series, and Maclaurin series expansions of functions.) Please note that it's been a while since I did any of this so I couldn't tell you, for example, exactly how to calculate the number of terms in the series you need to include guarantee an error that small enough to be negligible in a double-precision calculation.

For example, the Taylor/Maclaurin series expansion for e^x is this:

      +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

If you take all of the terms (k from 0 to infinity), this expansion is exact and complete (no error).

However, if you don't take all the terms going to infinity, but you stop after say 5 terms or 50 terms or whatever, you produce an approximate result that differs from the actual e^x function value by a remainder which is fairly easy to calculate.

The good news for exponentials is that it converges nicely and the terms of its polynomial expansion are fairly easy to code iteratively, so you might (repeat, MIGHT - remember, it's been a while) not even need to pre-calculate how many terms you need to guarantee your error is less than precision because you can test the size of the contribution at each iteration and stop when it becomes close enough to zero. In practice, I do not know if this strategy is viable or not - I'd have to try it. There are important details I have long since forgotten about. Stuff like: machine precision, machine error and rounding error, etc.

Also, please note that if you are not using e^x, but you are doing growth/decay with another base like 2^x or 10^x, the approximating polynomial function changes.

回复收藏 0 原文

别在捏我脸啦 2024-07-13 18:47:59

对于整数指数，将 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

The usual approach, to raise a to the b, for an integer exponent, goes something like this:

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

It is generally logarithmic in the size of the exponent. The algorithm is based on the invariant "a^b*result = a0^b0", where a0 and b0 are the initial values of a and b.

For negative or non-integer exponents, logarithms and approximations and numerical analysis are needed. The running time will depend on the algorithm used and what precision the library is tuned for.

Edit: Since there seems to be some interest, here's a version without the extra multiplication.

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

回复收藏 0 原文

梅窗月明清似水 2024-07-13 18:47:59

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

回复收藏 0 原文

绿萝 2024-07-13 18:47:59

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

回复收藏 0 原文

今天小雨转甜 2024-07-13 18:47:59

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

x * log2(e) = log2(1 + fraction) + exponent,     0 <= log2(1 + fraction) < 1

exponent = floor(x * log2(e))

1 + fraction = 2^(x * log2(e) - exponent) = e^((x * log2(e) - exponent) * ln2) = e^(x - exponent * ln2), 0 <= x - exponent * ln2 < ln2

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

x * log2(e) = log2(1 + fraction) + exponent,     0 <= log2(1 + fraction) < 1

exponent = floor(x * log2(e))

1 + fraction = 2^(x * log2(e) - exponent) = e^((x * log2(e) - exponent) * ln2) = e^(x - exponent * ln2), 0 <= x - exponent * ln2 < ln2

回复收藏 0 原文

~没有更多了~