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.
Best Answer
pow
works withdouble
numbers. These represent numbers of the form s * 2^e where s is a 53 bit integer. Thereforedouble
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 ify
can be large, but I assume it's always less than 64, making this issue moot.