How is pow() calculated in C?

Question

Our professor said that you can't calculate a^b if a<0 using pow() because pow() uses natural logarithms to calculate it (a^b=e^{b 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?

How precisely `pow` is calculated is a longer story, but the reason that it cannot be calculated when the base is negative and the exponent is non-integer is that the result would be a complex number that cannot be represented by the `double` that `pow` returns. There's `cpow` to handle that, though. — Wintermute, Nov 27 '16 at 00:45
Depends on your implementation, i.e. which compiler and C library you're using. — melpomene, Nov 27 '16 at 00:46
C99 7.12.7.4 (2) states: "The pow functions compute x raised to the power y. A domain error occurs if x is finite and negative and y is finite and not an integer value," so I don't think it depends on the implementation. `a` can be negative as long as `b` is infinite or integer. — Wintermute, Nov 27 '16 at 00:48
" have searched through math.h and further files, but was unable to find how the function is defined" - headers don't have definitions (`inline` functions are a specifial case). But I strongly doubt you could not find an implementation! There is at least one FOSS C library available. Is a search term for google like "c pow implementation" really **that** uncommon?? — too honest for this site, Nov 27 '16 at 00:52
With not optimized compiling, the result would be what inside the actual C library for your OS. With flags like `-O3`, the result is probably an optimized SSE/AVX implementation, which can be investigated by disassembling binary. — Kh40tiK, Nov 27 '16 at 00:53
Thanks @Wintermute . Now when you said that it's a complex number I'm even more curious to look how it works. — Szymon, Nov 27 '16 at 01:09
[Exponentiation is ambiguous](http://math.stackexchange.com/questions/317528/how-do-you-compute-negative-numbers-to-fractional-powers/317546#317546). The form of exponentiation that `double pow(double, double)` most naturally corresponds to is actually undefined for negative `a`, no matter what `b` is. — , Nov 27 '16 at 01:13
@Hurkyl: Are you saying that `(-5)^2` is ambiguous? It's really not. It's 25, no matter what way you look at it. — Dietrich Epp, Nov 27 '16 at 01:41
@Hurkyl What you have at start `(root c of (a^b)) = (root c of (a))^b)` of linked answer is true only under certain conditions. — Szymon, Nov 27 '16 at 02:04
@DietrichEpp: It's undefined in at least one way of looking at it. — , Nov 27 '16 at 02:15
@Hurkyl If you go `(-5)^2=(-5)^(4/2)=sqrt(-5)^4` than it's not directly computable. If you do such conversion you have to meet some conditions. If you calculate limits and get 0/0 you don't say you can't calculate the limit, you do conversions and get to equation for which you can easily calculate limit. With everything it depends on what you get at start. — Szymon, Nov 27 '16 at 02:28
@Szymon: I gave the conditions in my post where you can do that sort of arithmetic with discrete real exponentiation (your calculation is not one of them, since you have an even denominator). And for continuous exponentiation, the thing you start with, `(-5)^2` , is simply not in the domain. It may not seem important to you, but there really are contexts where continuity is more important than other concerns. Floating-point arithmetic, incidentally, is near the top of the list. I bet the C standard committee would not have allowed `pow(-5.0,2.0)` if C had `pow(double, int)` variation. — , Nov 27 '16 at 02:48
... but I wasn't trying to start a religious war here; the point of my comment was to increase the OP's awareness that there isn't a "one true answer" here; there really are different conventions (and to a lesser extent to point out there really are reasons for them too!). And yes, I'm aware the "kitchen sink" approach of throwing everything into the domain that you can reasonably (or even just semi-reasonably) justify is yet another conventions that some people like too, even if I personally don't like it. — , Nov 27 '16 at 02:58
@Hurkyl: There is, actually, one true answer for evaluating `pow(-5.0, 2.0)`. It will always be 25. You've cited an answer which talks about the differing mathematical definitions, but this is C with IEE-754 arithmetic, and the standards are quite clear on what the answer actually is. Yes, it's important to make people aware that different contexts have different conventions (for example, `log` is base-e, not 10) but in this case, we know which conventions are used from context. — Dietrich Epp, Nov 27 '16 at 03:22
…Another convention which appeared in this thread is `1.0/0.0` being infinity, which is quite clear from the IEEE-754 standard. Mathematically, it makes no sense, but we are not working with real numbers to begin with, so the rules are a little different. — Dietrich Epp, Nov 27 '16 at 03:24
@Szymon Maybe this [answer](http://stackoverflow.com/a/40519989/780717) can provide some insights as well. — njuffa, Nov 27 '16 at 18:03

Dietrich Epp · Accepted Answer · 2016-11-27T01:36:30.213

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.

Thanks. I guess I'll spend some time trying to understand it. On off-topic in e_pow.c there is this line `return y < 0 ? 1.0/0.0 : 0.0;` and got me thinking why return `1.0/0.0`? — Szymon, Nov 27 '16 at 01:47

score 6 · Answer 2 · answered Nov 27 '16 at 01:21

The second question (why does it return a domain error) is already covered in the comments, but adding for completeness: pow takes two real numbers and returns a real number. Applying a rational exponent on a negative number takes you out of the domain of real numbers into the domain of complex numbers, which the result of this function (a double) can't represent.

If you're curious about the actual implementation, well, there are many and it depends on many factors, such as architecture and level of optimisation. It's quite difficult to find one that reads easily, but 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.
 */

So, in short, the mechanism is as you described it and relies on calculating the logarithm first, but with many special cases that need to be accounted for.

For an enumeration of the special cases applicable to `pow()`, see also the ISO-C99 standard section F.9.4.4 (or the equivalent section of the latest ISO-C11 standard). — njuffa, Dec 06 '16 at 20:39

score 4 · Answer 3 · answered Nov 27 '16 at 01:20

Assuming an x86 series processor, pow is the equivalent of

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

Where exp2 and log2 are CPU primitives for the exponential and logarithm operations in base 2.

Different CPUs inherently have different implementations.

In theory if you didn't have pow you could write:

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

but this loses precision over the native version due to accumulative roundoff.

And Dietrich Epp revealed I missed a bunch of special cases. Nevertheless I have something to say about roundoff that should be allowed to stand.

Fails for half of all `double`: when `base <= 0` - not really a _special_ case. — chux - Reinstate Monica, Jan 05 '19 at 00:49

score 2 · Answer 4 · answered Nov 27 '16 at 02:06

pow does work 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 involves the y-th root of x. For instance, when you try to calculate a^1/2 you are actually looking for the square root of a.

So what happens if you have a negative base and a non-integer exponent? You get a y-th root of a negative number, which yields is a complex non-real number. pow() does not work with complex numbers, so it will probably return NaN.

How is pow() calculated in C?

4 Answers4

The C Standard

Linked