How is pow () calculated in C?

Question

Our professor said that you cannot calculate ^b if a <0 uses pow() , because pow() uses the natural logarithms to calculate it (a ^b = e ^{b ln a} ), and since it is undefined for negative numbers, it does not can be calculated. I tried and it works as long as b is an integer.

I looked at math.h and other files, but could not find how the function is defined and what it uses to calculate. I also tried searching the Internet, but without any success. The Stack Overflow section has similar questions: here and here (for C #). (the latter is good, but I could not find the source code.)

So the question is, how are pow() actually computed in C? And why does it return a domain error when the base is finite and negative, and the exponent is finite and not an integral?

Answer 1

If you are interested in how the pow function can be implemented in practice, you can look at the source code. There is a kind of “ability” to search through unfamiliar (and large) codebases to find the section you are looking for, and it would be nice to get some practice.

One implementation of the C library is glibc, which has mirrors on GitHub. I did not find an official mirror, but an unofficial mirror is located at https://github.com/lattera/glibc

First, consider the math/w_pow.c , which has a promising name. It contains the __pow function, which calls __ieee754_pow , which we can find in sysdeps/ieee754/dbl-64/e_pow.c (remember that not all IEEE-754 systems, so it makes sense that the IEEE-754 math code is in its own catalog).

It begins with several special cases:

 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; 

A little further you will find a comment thread

 /* if x<0 */ 

What brings us to

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

So you can see that for negative x and integer y version of glibc pow will calculate pow(-x,y) and then make the result negative if y is odd.

This is not the only way to do something, but I assume this is common to many implementations. You can see that pow filled with special cases. This is often found in library math functions, which should work correctly with unfriendly inputs such as denormals and infinity.

The pow function is especially difficult to read because it is highly optimized code that performs bit-twisting floating point numbers.

Standard c

The C standard (n1548 §7.12.7.4) refers to pow :

A domain error occurs if x is finite and negative, and y is finite, not an integer value.

So, according to the C standard, negative x should work.

There is also Appendix F, which gives much more stringent restrictions on how pow works on IEEE-754 / IEC-60559 systems.

Answer 2

The second question (why it returns a domain error) has already been discussed in the comments, but the addition is for completeness: pow takes two real numbers and returns a real number. Using a rational indicator on a negative number takes you from the real number region to the complex number region, the result of which cannot be represented by this function (double).

If you are interested about a real implementation, then there are many of them, and this depends on many factors, such as architecture and level of optimization. It's hard to find one that reads easily, but FDLIBM (Freely Distributable LIBM) has one that has at least 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. */ 

So, in short, the mechanism, as you described it, first calculates the logarithm, but with many special cases that need to be considered.

Answer 3

Assuming an x86 series processor, pow is the equivalent

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

Where exp2 and exp2 are CPU primitives for exponential and logarithmic operations in base 2.

Different processors inherently have different implementations.

In theory, if you did not have pow , you could write:

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

but it loses accuracy over the original version due to cumulative rounding.

And Dietrich Epp testified that I missed a bunch of special occasions. However, I need to say something about rounding, which should be allowed to stand.

Answer 4

pow works for negative numbers. It just doesn't work when the base is negative, and the exponent is not an integer.

A number in the form a ^{x / y} actually includes the yth root of x. For example, when you try to calculate ^1/2 you are really looking for the square root of.

So what happens if you have a negative base and a non-integer indicator? You get the yth root of a negative number, which gives a complex unreal number. pow() does not work with complex numbers, so it will probably return NaN.