我们的教授说你不能计算一个b,如果a <0使用pow()因为pow()使用自然对数来计算它(a b = e b ln a),并且由于它未定义为负数,所以无法计算.我尝试了它,只要b是一个整数就可以工作.
我已经搜索了math.h更多文件,但无法找到函数的定义方式以及它用于计算的内容.我也试过在互联网上搜索,但没有任何成功.Stack Overflow就在这里和这里有类似的问题(对于C#).(最后一个很好,但我找不到源代码.)
那么问题是如何pow()在C中实际计算？当基数是有限的且负数且指数是有限的和非整数时,为什么它会返回域错误？
Die*_*Epp 13
如果您对如何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标准(n1548§7.12.7.4)有这样的说法pow:
如果x是有限的并且是负的且y是有限的而不是整数值,则会发生域错误.
因此,根据C标准,否定x 应该起作用.
还有附录F的问题,它对如何pow在IEEE-754/IEC-60559系统上工作提出了更严格的限制.
第二个问题(为什么它会返回域错误)已在评论中介绍,但添加完整性:pow需要两个实数并返回一个实数.在负数上应用有理指数会使您从实数域进入复数域,这个函数的结果(双精度)无法表示.
如果您对实际实现感到好奇,那么有很多,这取决于许多因素,例如架构和优化级别.找到一个易于阅读的内容非常困难,但FDLIBM(Freely Distributable LIBM)在评论中至少有一个很好的解释:
/* __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.
*/
因此,简而言之,该机制与您所描述的一样,并且首先依赖于计算对数,但需要考虑许多特殊情况.
假设 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));
}
但由于累积舍入，这会失去本机版本的精度。
迪特里希·埃普透露我错过了很多特殊情况。不过，关于舍入，我有一些话要说，应该允许保留。