c - 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
", pow(17, 12));
    printf("17^13 = %lf
", pow(17, 13));
    printf("17^14 = %lf
", 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.

question from:https://stackoverflow.com/questions/21452711/pow-seems-to-be-out-by-one-here

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.