在 C 中如何计算 pow()？答案

【问题标题】：How is pow() calculated in C?在 C 中如何计算 pow()？
【发布时间】：2017-04-11 00:24:21
【问题描述】：

我们的教授说如果apow() 无法计算a^b，因为pow() 使用自然对数来计算它（a^b= e^{b ln a}) 并且由于它对于负数未定义，因此无法计算。我试过了，只要 b 是整数，它就可以工作。

我搜索了math.h 和其他文件，但无法找到该函数的定义方式以及它用于计算的内容。我也尝试搜索互联网，但没有任何成功。在 Stack Overflow 上也有类似的问题 here 和 here（对于 C#）。（最后一个不错，就是找不到源代码。）

所以问题是pow() 在 C 中实际上是如何计算的？为什么当基数是有限且负数而指数是有限且非整数时，它会返回域错误？

【问题讨论】：

pow 的计算精度是一个较长的故事，但是当底数为负数且指数为非整数时无法计算它的原因是结果将是一个复数不能由 pow 返回的 double 表示。不过，有 cpow 来处理这个问题。
取决于您的实现，即您使用的编译器和 C 库。
C99 7.12.7.4 (2) 状态：“pow 函数计算 x 的 y 次幂。如果 x 是有限的负数并且 y 是有限的而不是整数值，则会发生域错误，”所以我认为这不取决于实施。 a 可以是负数，只要 b 是无限的或整数。
" 已搜索 math.h 和其他文件，但无法找到函数的定义方式" - 标头没有定义（inline 函数是特殊情况）。但我强烈怀疑你找不到实现！至少有一个 FOSS C 库可用。像“c pow implementation”这样的谷歌搜索词真的那不常见吗？？
如果没有优化编译，结果将是您操作系统的实际 C 库中的内容。使用 -O3 这样的标志，结果可能是优化的 SSE/AVX 实现，可以通过反汇编二进制来研究。

标签： c pow math.h

【解决方案1】：

如果您对pow 函数在实践中的实现方式感到好奇，您可以查看源代码。有一种“诀窍”可以在不熟悉的（和大型的）代码库中搜索以找到您要查找的部分，并且最好进行一些练习。

C 库的一个实现是 glibc，它在 GitHub 上有镜像。没找到官方镜像，非官方镜像在https://github.com/lattera/glibc

我们首先看一下math/w_pow.c 文件，它的名字很有前途。它包含一个函数__pow，它调用__ieee754_pow，我们可以在sysdeps/ieee754/dbl-64/e_pow.c 中找到它（请记住，并非所有系统都是IEEE-754，因此IEEE-754 数学代码位于其自己的目录中是有意义的）。

从一些特殊情况开始：

if (y == 1.0) return x;
if (y == 2.0) return x*x;
if (y == -1.0) return 1.0/x;
if (y == 0) return 1.0;

再往下一点你会发现一个带有评论的分支

/* if x<0 */

这导致我们去

return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */

因此您可以看到，对于负数 x 和整数 y，glibc 版本的 pow 将计算 pow(-x,y)，然后如果 y 是奇数则将结果设为负数。

这不是唯一的做事方式，但我猜这对许多实现来说很常见。您可以看到pow 充满了特殊情况。这在库数学函数中很常见，它们应该与非正规和无穷大等不友好的输入正常工作。

pow 函数特别难以阅读，因为它是高度优化的代码，对浮点数进行位旋转。

C 标准

C 标准（n1548 §7.12.7.4）对pow 有这样的说法：

如果 x 是有限的负数并且 y 是有限的而不是整数值，则会发生域错误。

因此，根据 C 标准，否定 x 应该起作用。

还有附录 F 的问题，它对 pow 如何在 IEEE-754 / IEC-60559 系统上工作给出了更严格的限制。

【讨论】：

谢谢。我想我会花一些时间试图理解它。在 e_pow.c 的题外话上有这行 return y < 0 ? 1.0/0.0 : 0.0; 让我思考为什么要返回 1.0/0.0？
1.0/0.0 是无穷大。

【解决方案2】：

第二个问题（为什么它返回域错误）已经包含在 cmets 中，但为了完整性添加：pow 接受两个实数并返回一个实数。对负数应用有理指数会将您从实数域带入复数域，此函数的结果 (a double) 无法表示。

如果您对实际实现感到好奇，那么有很多，这取决于很多因素，例如架构和优化级别。很难找到一个容易阅读的，但是 FDLIBM (Freely Distributable LIBM) has one which has at least has a good explanation in the comments:

/* __ieee754_pow(x,y) return x**y
 *
 *            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)
 *
 * Special cases:
 *  1.  (anything) ** 0  is 1
 *  2.  (anything) ** 1  is itself
 *  3.  (anything) ** NAN is NAN
 *  4.  NAN ** (anything except 0) is NAN
 *  5.  +-(|x| > 1) **  +INF is +INF
 *  6.  +-(|x| > 1) **  -INF is +0
 *  7.  +-(|x| < 1) **  +INF is +0
 *  8.  +-(|x| < 1) **  -INF is +INF
 *  9.  +-1         ** +-INF is NAN
 *  10. +0 ** (+anything except 0, NAN)               is +0
 *  11. -0 ** (+anything except 0, NAN, odd integer)  is +0
 *  12. +0 ** (-anything except 0, NAN)               is +INF
 *  13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
 *  14. -0 ** (odd integer) = -( +0 ** (odd integer) )
 *  15. +INF ** (+anything except 0,NAN) is +INF
 *  16. +INF ** (-anything except 0,NAN) is +0
 *  17. -INF ** (anything)  = -0 ** (-anything)
 *  18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
 *  19. (-anything except 0 and inf) ** (non-integer) is NAN
 *
 * Accuracy:
 *  pow(x,y) returns x**y nearly rounded. In particular
 *          pow(integer,integer)
 *  always returns the correct integer provided it is 
 *  representable.
 *
 * Constants :
 * The hexadecimal values are the intended ones for the following 
 * constants. The decimal values may be used, provided that the 
 * compiler will convert from decimal to binary accurately enough 
 * to produce the hexadecimal values shown.
 */

因此，简而言之，该机制正如您所描述的那样，并且首先依赖于计算对数，但需要考虑许多特殊情况。

【讨论】：

有关适用于pow() 的特殊情况的枚举，另请参见 ISO-C99 标准部分 F.9.4.4（或最新 ISO-C11 标准的等效部分）。跨度>

【解决方案3】：

假设 x86 系列处理器，pow 相当于

double pow(double base, double exp)
{
   return exp2(exp * log2(base));
}

其中 exp2 和 log2 是用于以 2 为底的指数和对数运算的 CPU 原语。

不同的 CPU 天生就有不同的实现。

理论上，如果你没有pow，你可以写：

double pow(double base, double exponent)
{
   return exp(exponent * log(base));
}

但是由于累积舍入，这会比原生版本失去精度。

Dietrich Epp 透露我错过了一些特殊情况。不过，我有话要说，关于舍入应该被允许。

【讨论】：

一半的double 失败：当base <= 0 - 不是真正的特殊案例。

【解决方案4】：

pow 对负数有效。当底数为负且指数不是整数时，它就不起作用了。

a^x/y 形式的数实际上涉及 x 的第 y 根。例如，当您尝试计算 a^1/2 时，实际上是在寻找 a 的平方根。

如果你有一个负数和一个非整数指数会发生什么？你得到一个负数的第 y 根，它产生一个复数的非实数。 pow() 不适用于复数，因此它可能会返回 NaN。

【讨论】：