CODEWITHSUNDEEP
cgccprecisionpow
igorludi
igorludi

Reputation: 1596

pow numeric error in c

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.

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;
}

Upvotes: 3

Views: 203

Answers (1)

Bathsheba
Bathsheba

Reputation: 234875

(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.

Upvotes: 2

Related Questions