Search code examples
cpowmath.h

How is pow() calculated in C?

Our professor said that you can't calculate ab if a<0 using pow() because pow() uses natural logarithms to calculate it (ab=eb ln a) and since it's undefined for negative numbers it can't be calculated. I tried it and it works as long as b is an integer.

I have searched through math.h and further files, but was unable to find how the function is defined and what it uses to calculate. I also tried searching the internet, but without any success. There are similar questions on Stack Overflow right here and here (for C#). (the last one is good, but I was unable to find sourcecode.)

So the question is how is pow() actually calculated in C? And why does it return a domain error when the base is finite and negative and the exponent is finite and non-integral?

Solution

  • If you're curious how the pow function might be implemented in practice, you can look at the source code. There is a kind of "knack" to searching through unfamiliar (and large) codebases to find the section you are looking for, and it's good to get some practice.

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

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

    It starts with a few 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 farther down you find a branch with a comment

    /* if x<0 */

    Which leads us to

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

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

    This is not the only way to do things, but my guess is that this is common to many implementations. You can see that pow is full of special cases. This is common in library math functions, which are supposed to work correctly with unfriendly inputs like denormals and infinity.

    The pow function is especially hard to read because it is heavily-optimized code which does bit-twiddling on floating-point numbers.

    The C Standard

    The C standard (n1548 §7.12.7.4) has this to say about pow:

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

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

    There is also the matter of appendix F, which gives much tighter constraints on how pow works on IEEE-754 / IEC-60559 systems.