pow() seems to be out by one here
c, floating-point, floating-point-precision, pow
Solution
`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.
Problem
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.