Szymon - 1 year ago 61
C Question

# 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?

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 */
``````

``````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.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download