Will (int)pow(n,m) be wrong for some positive integers n,m?

Question

Assuming n and m are positive integers, and n^m is within the range of an integer, will (int)pow(n,m) ever give a wrong answer?

I have tried many n for m=2 and have not gotten any wrong answers so far.

Answer 1

The C standard does not impose any requirements on the accuracy of floating point arithmetic. The accuracy is implementation-defined which means that implementations are required to document it. However, implementations are left with a significant "out": (§5.2.4.2.2 paragraph 6, emphasis added.)

The accuracy of the floating-point operations (+, -, *, /) and of the library functions in <math.h> and <complex.h> that return floating-point results is implementation-defined, as is the accuracy of the conversion between floating-point internal representations and string representations performed by the library functions in <stdio.h>, <stdlib.h>, and <wchar.h>. The implementation may state that the accuracy is unknown.

And, indeed, gcc takes advantage of this by specifying that the accuracy is unknown. Nonetheless, the accuracy of glibc computations is pretty good, even if not guaranteed.

The MS libc implementation is known to occasionally produce an error of 1ULP for the pow function with integer arguments, resulting in an incorrect value if the result of the pow operation is simply truncated to an int. (I couldn't find any specification in the Visual Studio documentation about floating point accuracy, but I believe that the list of SO questions below provides evidence of my assertion.)

On x86 architectures, most implementations make some attempt to implement IEEE 754, since the native floating point representation conforms. However, until the 2008 revision, IEEE-754 only required correctly-rounded results from +, -, *, / and sqrt. Since the revision, it recommends that a number of other functions return correctly-rounded results, but all of these recommendations are optional, and few math libraries implement all of them.

If you really want to use pow to compute integer powers of integers, it is advisable (and easy) to use lround(pow(n, m)) instead of (long)(pow(n, m)), which will round the result to the nearest integer, rather than optimistically relying on the error being positive. That should give the correct integer value for results up to 2⁵² (with IEEE-754 doubles) if pow is within 1ULP. Between 2⁵² and 2⁵³, a 1ULP error would be 0.5, which will sometimes round to the wrong integer. Beyond 2⁵³ not all integers are representable as doubles.

SO is actually full of questions resulting from this particular problem. See:

• Reversing a five digit number with POW function in C
• C's pow() function as per header file <math.h> not working properly
• Strange behaviour of the pow function
• Why am I getting unexpected output when using floor with pow?
• Why pow(10,5) = 9,999 in C++
• return value of pow() gets rounded down if assigned to an integer (found by Lưu Vĩnh Phúc)
• wrong output by power function - C (also)

and undoubtedly many more.

Answer 2

Some implementations evaluate pow(x,y) as exp(y*log(x)) in any case, and others use a square-multiply powering for integral exponents.

glibc for example has 'C' implementations, as well as platform-specific implementations with lookup tables, polynomial approximations, etc. Many Sun/BSD derivatives appear to handle integral exponents as a special case.

More importantly, IEEE-754 doesn't require it. Nor does it specify the required accuracy of the result. glibc documents its maximum ulp errors, though this page might not be up to date.

In summary, don't assume an exact integral result, even if pow(n,m) has an exact floating-point representation.