Approximation of pow() in Java, C and C++

Author: Martin Ankerl, 30/07/08 Keywords: Optimization, Java, C, CPP Abstract: For many application of the pow() method, an approximate value suffice. This text shows how to gain a 41 times speed increase by using an approximated value that is much easier to calculate. It is a simple "how to" for the languages Java, C and C++.

subscribe to my RSS feed

Approximation of pow() in Java, C and C++

Approximation of pow() in Java

public static double pow(final double a, final double b) {
    final int x = (int) (Double.doubleToLongBits(a) >> 32);
    final int y = (int) (b * (x - 1072632447) + 1072632447);
    return Double.longBitsToDouble(((long) y) << 32);
}

This is really very compact. The calculation only requires 2 shifts, 1 mul, 2 add, and 2 register operations. That's it! In my tests it usually within an error margin of 5% to 12%, in extreme cases sometimes up to 25%. A careful analysis is left as an exercise for the reader. This is very usable for in e.g. metaheuristics or neural nets.

I use Linux, Java 1.6.0-b105 with the server VM, and execute the benchmark with this command: sudo nice -n -20 java -cp . -server PowTest

The approximation is 27 times faster than Math.pow() on my Pentium-M. On a Pentium 4 it is 41 times faster. Unfortunately, micro benchmarks are difficult to do in Java, so your mileage may vary.

Approximation of pow() in C and C++

double pow(double a, double b) {
    int tmp = (*(1 + (int *) &a));
    int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447);
    double p = 0.0;
    *(1 + (int * ) &p) = tmp2;
    return p;
}

This version is 7.8 times faster than pow() from the standard library.

WARNING! you HAVE to use the -fno-strict-aliasing option, or this does not work!

How the Approximation was Developed

Approximation of e^x

public static double exp(double val) {
    final long tmp = (long) (1512775 * val + (1072693248 - 60801));
    return Double.longBitsToDouble(tmp << 32);
}

Use Exponential Functions for a^b

1. Take exponential a^b = e^(ln(a^b))
2. Extract b a^b = e^(ln(a)*b)

Now we have expressed the pow calculation with e^x and ln(x). We already have the e^x approximation, but no good ln(x).

Approximation of ln(x)

That's not too difficult. Have a look at the above code, we now take the output and move backwards to undo the calculation. First reverse the shift:

final double tmp = (Double.doubleToLongBits(val) >> 32);

tmp = (1512775 * val + (1072693248 - 60801))

for val:

1. The original formula tmp = (1512775 * val + (1072693248 - 60801))
2. Perform subtraction tmp = 1512775 * val + 1072632447
3. Bring value to other side tmp - 1072632447 = 1512775 * val
4. Divide by factor (tmp - 1072632447) / 1512775 = val
5. Finally, val on the left side val = (tmp - 1072632447) / 1512775

Voíla, now we have a nice approximation of ln(x):

public double ln(double val) {
    final double x = (Double.doubleToLongBits(val) >> 32);
    return (x - 1072632447) / 1512775;
}

Combine Both Approximations

public static double pow1(final double a, final double b) {
    // calculate ln(a)
    final double x = (Double.doubleToLongBits(a) >> 32);
    final double ln_a = (x - 1072632447) / 1512775;

    // ln(a) * b
    final double tmp1 = ln_a * b;

    // e^(ln(a) * b)
    final long tmp2 = (long) (1512775 * tmp1 + (1072693248 - 60801));
    return Double.longBitsToDouble(tmp2 << 32);
}

Between the two shifts, we can simply insert the tmp1 calculation into the tmp2 calculation to get

public static double pow2(final double a, final double b) {
    final double x = (Double.doubleToLongBits(a) >> 32);
    final long tmp2 = (long) (1512775 * (x - 1072632447) / 1512775 * b + (1072693248 - 60801));
    return Double.longBitsToDouble(tmp2 << 32);
}

Now simplify tmp2 calculation:

1. The original formula tmp2 = (1512775 * (x - 1072632447) / 1512775 * b + (1072693248 - 60801))
2. We can drop the factor 1512775 tmp2 = (x - 1072632447) * b + (1072693248 - 60801)
3. And finally, calculate the substraction tmp2 = b * (x - 1072632447) + 1072632447

The Result

public static double pow(final double a, final double b) {
    final int tmp = (int) (Double.doubleToLongBits(a) >> 32);
    final int tmp2 = (int) (b * (tmp - 1072632447) + 1072632447);
    return Double.longBitsToDouble(((long) tmp2) << 32);
}

This concludes my little tutorial on micro optimization of the pow() function. If you have come this far, I congratulate your persistence :-)

Addendum

a between 0 and 1000, b between 0 and 5:
Worst case:
a=512.0125338006894
b=4.914054794454942
a^b approximation: 2.5571362865152E13
a^b Math.pow(): 2.0585114388503066E13
Errors is 19.499345822682237 %

Average error is 4.021374964371438 %

---
a between 0 and 100, b between 0 and 3:
Worst case:
a=64.00103767757574
b=2.8915318496742626
a^b approximation: 191223.125
a^b Math.pow(): 166973.39656298532
Error is 12.681378592162784 %

Average error is 2.7778168699408558 %

I am using this (a modified form) in an Ant Colony Optimization algorithm. There the core algorithms needs to calculate a^b * c^d very often. I have done comparison of the optimization quality when using this optimization and when using Math.pow(), the error does not have any measurable negative effect in the examples I am using. And the optimization got a hell of a lot faster, so even if there is a slight decrease in optimization quality this is more than compensated by the higher number of iterations I am able to do now.

Comments

• Artur Biesiadowski
  Unfortunately, you equation fails with HUGE errors when powers go up. You can find quite good middle ground by combining two approaches
  • compute integer part of power by just multiplying the number by itself n times
  • then use your equation for rest of the power (0-1) and multiply both results.
  With that approach, it is bit slower, but max error seems to be bound to around 5% regardless of the power. Maybe there is another bit-related trick to get good approximation of a^b where b is integer, which could be combined together to have constant time?
• Amphi
  Another thing one might consider are look up tables. E.g. if you only need 256 different pow results, a LUT might be the better choice. (Where "better" means only a few percent faster than an approximation but far easier to understand.)
  E.g. if you're working with colors you've typically that 0-255 range and one of the parameters is a fixed value.
• Kasper Graversen
  Good suggestion. You touch upon the ever lasting problem in computer science– the trade-off between memory and speed.

Help spread the word