3

I'm wondering where does the numeric error happen, in what layer. Let me explain using an example:

int p = pow(5, 3);
printf("%d", p);

I've tested this code on various HW and compilers (VS and GCC) and some of them print out 124, and some 125.

  • On the same HW (OS) i get different results in different compilers (VS and GCC).
  • On the different HW(OS) I get different results in the same compiler (cc (GCC) 4.8.1).

AFAIK, pow computes to 124.99999999 and that gets truncated to int, but where does this error happen? Or, in other words, where does the correction happen (124.99->125)

Is it a compiler-HW interaction?

//****** edited:

Here's an additional snippet to play with (keep an eye on p=5, p=18, ...):

#include <stdio.h>
#include <math.h>

int main(void) {
    int p;
    for (p = 1; p < 20; p++) {
        printf("\n%d %d %f %f", (int) pow(p, 3), (int) exp(3 * log(p)), pow(p, 3), exp(3 * log(p)));
    }
    return 0;
}
igorludi
  • 1,519
  • 2
  • 18
  • 31
  • 2
    Floating point values on computers are a complex matter, and results of operations on them depend on many things. For example, the hardware support (if any), rounding algorithms used by the different software (compilers, standard libraries) or hardware, and many many other things. – Some programmer dude Oct 19 '15 at 08:48
  • 1
    Also be aware then when using constant expressions the compiler often will [use builtin functions](http://stackoverflow.com/a/24294632/1708801) and so you may also want to try this using `-fno-builtin` to see if that effects the results. – Shafik Yaghmour Oct 19 '15 at 09:21
  • I believe it's not, I've just added a snippet using a for loop and the same thing happens. – igorludi Oct 19 '15 at 09:27
  • 2
    Why are truncating to `int`? Why are you using `pow` with an integral exponent? These are the real questions to be asked. – David Heffernan Oct 19 '15 at 09:32
  • True, but as I said - I'm aware this is not the best practice and I'm surely not advocating this. As a side note, the question popped up when one student tried to compute whether a number is an Armstrong number (http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap04/arms.html). It is an honest mistake, especially for a student. – igorludi Oct 19 '15 at 10:02

1 Answers1

2

(First note that for an IEEE754 double precision floating point type, all integers up to the 53rd power of 2 can be represented exactly. Blaming floating point precision for integral pow inaccuracies is normally incorrect).

pow(x, y) is normally implemented in C as exp(y * log(x)). Hence it can "go off" for even quite small integral cases.

For small integral cases, I normally write the computation long-hand, and for other integral arguments I use a 3rd party library. Although a do-it-yourself solution using a for loop is tempting, there are effective optimisations that can be done for integral powers that such a solution might not exploit.

As for the observed different results, it could be down to some of the platforms using an 80 bit floating point intermediary. Perhaps some of the computations then are above 125 and others are below that.

Bathsheba
  • 231,907
  • 34
  • 361
  • 483
  • OK, I get that, and I'm aware that this isn't the best practice (pow->int). I also get that pow creates rounding errors, but I don't get where does this error occur, as I said in the question. – igorludi Oct 19 '15 at 09:01
  • 2
    As for those looking for effective int pow implementation: http://stackoverflow.com/questions/101439/the-most-efficient-way-to-implement-an-integer-based-power-function-powint-int – igorludi Oct 19 '15 at 09:02
  • 1
    `log(5)` is a `double` which **can't** be represented exactly. `exp(3 log (5))` will, most likely, **not** be a whole number. – Bathsheba Oct 19 '15 at 09:02
  • 1
    Try this code: #include #include int main(void) { int p; for (p = 1; p < 20; p++) { printf("\n%d %d %f %f", (int) pow(p, 3), (int) exp(3 * log(p)), pow(p, 3), exp(3 * log(p))); } return 0; } Using GCC, I get: 124 124 125.000000 125.000000 Using VS, I get: 125 124 125.000000 125.000000 I'm not sure that pow is implemented that way, or VS performs corrections on pow. Either way, what still puzzles me it that I'd get different results for GCC on another machine. There are various forces at play here... – igorludi Oct 19 '15 at 09:12
  • I suspect the differences are down to 80 bit intermediaries. And perhaps the C standard library is implementing `pow` in a different way on platforms. There are various compiler settings you can tweak to change the former. – Bathsheba Oct 19 '15 at 09:14
  • As I pointed out in my comment above, since the example using constant expressions bultins may be used as well – Shafik Yaghmour Oct 19 '15 at 09:24