Numerical accuracy of pow(a/b,x) vs pow(b/a,-x)

Question

Is there a difference in accuracy between pow(a/b,x) and pow(b/a,-x)? If there is, does raising a number less than 1 to a positive power or a number greater than 1 to a negative power produce more accurate result?

Edit: Let's assume x86_64 processor and gcc compiler.

Edit: I tried comparing using some random numbers. For example:

printf("%.20f",pow(8.72138221/1.761329479,-1.51231)) // 0.08898783049228660424
printf("%.20f",pow(1.761329479/8.72138221, 1.51231)) // 0.08898783049228659037

So, it looks like there is a difference (albeit minuscule in this case), but maybe someone who knows about the algorithm implementation could comment on what the maximum difference is, and under what conditions.

Answer 1

Here's one way to answer such questions, to see how floating-point behaves. This is not a 100% correct way to analyze such question, but it gives a general idea.

Let's generate random numbers. Calculate v0=pow(a/b, n) and v1=pow(b/a, -n) in float precision. And calculate ref=pow(a/b, n) in double precision, and round it to float. We use ref as a reference value (we suppose that double has much more precision than float, so we can trust that ref can be considered the best possible value. This is true for IEEE-754 for most of the time). Then sum the difference between v0-ref and v1-ref. The difference should calculated with "the number of floating point numbers between v and ref".

Note, that the results may be depend on the range of a, b and n (and on the random generator quality. If it's really bad, it may give a biased result). Here, I've used a=[0..1], b=[0..1] and n=[-2..2]. Furthermore, this answer supposes that the algorithm of float/double division/pow is the same kind, have the same characteristics.

For my computer, the summed differences are: 2604828 2603684, it means that there is no significant precision difference between the two.

Here's the code (note, this code supposes IEEE-754 arithmetic):

#include <cmath>
#include <stdio.h>
#include <string.h>

long long int diff(float a, float b) {
    unsigned int ai, bi;
    memcpy(&ai, &a, 4);
    memcpy(&bi, &b, 4);
    long long int diff = (long long int)ai - bi;
    if (diff<0) diff = -diff;
    return diff;
}

int main() {
    long long int e0 = 0;
    long long int e1 = 0;
    for (int i=0; i<10000000; i++) {
        float a = 1.0f*rand()/RAND_MAX;
        float b = 1.0f*rand()/RAND_MAX;
        float n = 4.0f*rand()/RAND_MAX - 2.0f;

        if (a==0||b==0) continue;

        float v0 = std::pow(a/b, n);
        float v1 = std::pow(b/a, -n);
        float ref = std::pow((double)a/b, n);

        e0 += diff(ref, v0);
        e1 += diff(ref, v1);
    }

    printf("%lld %lld\n", e0, e1);
}

Answer 2

... between pow(a/b,x) and pow(b/a,-x) ... does raising a number less than 1 to a positive power or a number greater than 1 to a negative power produce more accurate result?

Whichever division is more arcuate.

---

Consider z = x^y = 2^{y * log2(x)}.

Roughly: The error in y * log2(x) is magnified by the value of z to form the error in z. x^y is very sensitive to the error in x. The larger the |log2(x)|, the greater concern.

In OP's case, both pow(a/b,p) and pow(b/a,-p), in general, have the same y * log2(x) and same z and similar errors in z. It is a question of how x, y are formed:

---

a/b and b/a, in general, both have the same error of +/- 0.5*unit in the last place and so both approaches are of similar error.

Yet with select values of a/b vs. b/a, one quotient will be more exact and it is that approach with the lower pow() error.

pow(7777777/4,-p) can be expected to be more accurate than pow(4/7777777,p).

Lacking assurance about the error in the division, the general rule applies: no major difference.