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;
}
5
  • 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. Oct 19, 2015 at 8:48
  • 1
    Also be aware then when using constant expressions the compiler often will use builtin functions and so you may also want to try this using -fno-builtin to see if that effects the results. Oct 19, 2015 at 9:21
  • I believe it's not, I've just added a snippet using a for loop and the same thing happens.
    – igorludi
    Oct 19, 2015 at 9:27
  • 2
    Why are truncating to int? Why are you using pow with an integral exponent? These are the real questions to be asked. Oct 19, 2015 at 9: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 (cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap04/arms.html). It is an honest mistake, especially for a student.
    – igorludi
    Oct 19, 2015 at 10:02

1 Answer 1

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.

6
  • 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, 2015 at 9:01
  • 2
    As for those looking for effective int pow implementation: stackoverflow.com/questions/101439/…
    – igorludi
    Oct 19, 2015 at 9: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, 2015 at 9:02
  • 1
    Try this code: #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; } 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, 2015 at 9: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, 2015 at 9:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.