pow() seems to be out by one here

Question

What's going on here:

#include <stdio.h>
#include <math.h>
int main(void) {
    printf("17^12 = %lf\n", pow(17, 12));
    printf("17^13 = %lf\n", pow(17, 13));
    printf("17^14 = %lf\n", pow(17, 14));
}

I get this output:

17^12 = 582622237229761.000000
17^13 = 9904578032905936.000000
17^14 = 168377826559400928.000000

13 and 14 do not match with wolfram alpa cf:

12: 582622237229761.000000
    582622237229761

13: 9904578032905936.000000
    9904578032905937

14: 168377826559400928.000000
    168377826559400929

Moreover, it's not wrong by some strange fraction - it's wrong by exactly one!

If this is down to me reaching the limits of what pow() can do for me, is there an alternative that can calculate this? I need a function that can calculate x^y, where x^y is always less than ULLONG_MAX.

Answer 1

pow works with double numbers. These represent numbers of the form s * 2^e where s is a 53 bit integer. Therefore double can store all integers below 2^53, but only some integers above 2^53. In particular, it can only represent even numbers > 2^53, since for e > 0 the value is always a multiple of 2.

17^13 needs 54 bits to represent exactly, so e is set to 1 and hence the calculated value becomes even number. The correct value is odd, so it's not surprising it's off by one. Likewise, 17^14 takes 58 bits to represent. That it too is off by one is a lucky coincidence (as long as you don't apply too much number theory), it just happens to be one off from a multiple of 32, which is the granularity at which double numbers of that magnitude are rounded.

For exact integer exponentiation, you should use integers all the way. Write your own double-free exponentiation routine. Use exponentiation by squaring if y can be large, but I assume it's always less than 64, making this issue moot.

Answer 2

The numbers you get are too big to be represented with a double accurately. A double-precision floating-point number has essentially 53 significant binary digits and can represent all integers up to 2^53 or 9,007,199,254,740,992.

For higher numbers, the last digits get truncated and the result of your calculation is rounded to the next number that can be represented as a double. For 17^13, which is only slightly above the limit, this is the closest even number. For numbers greater than 2^54 this is the closest number that is divisible by four, and so on.

Answer 3

If your input arguments are non-negative integers, then you can implement your own pow.

Recursively:

unsigned long long pow(unsigned long long x,unsigned int y)
{
    if (y == 0)
        return 1;
    if (y == 1)
        return x;
    return pow(x,y/2)*pow(x,y-y/2);
}

Iteratively:

unsigned long long pow(unsigned long long x,unsigned int y)
{
    unsigned long long res = 1;
    while (y--)
        res *= x;
    return res;
}

Efficiently:

unsigned long long pow(unsigned long long x,unsigned int y)
{
    unsigned long long res = 1;
    while (y > 0)
    {
        if (y & 1)
            res *= x;
        y >>= 1;
        x *= x;
    }
    return res;
}

Answer 4

A small addition to other good answers: under x86 architecture there is usually available x87 80-bit extended format, which is supported by most C compilers via the long double type. This format allows to operate with integer numbers up to 2^64 without gaps.

There is analogue of pow() in <math.h> which is intended for operating with long double numbers - powl(). It should also be noticed that the format specifier for the long double values is other than for double ones - %Lf. So the correct program using the long double type looks like this:

#include <stdio.h>
#include <math.h>
int main(void) {
    printf("17^12 = %Lf\n", powl(17, 12));
    printf("17^13 = %Lf\n", powl(17, 13));
    printf("17^14 = %Lf\n", powl(17, 14));
}

As Stephen Canon noted in comments there is no guarantee that this program should give exact result.