10

I've read that cmath calculates pow(a,b) by performing exp(b*log(a)). This should not be used when b is an integer, since it slows down calculations a lot. What alternatives are there when

  1. calculating a lot of successive pow()s with the same constant a
  2. it is known beforehand that b will definitely be an integer?

I am looking for fast alternatives which are efficient in these particular scenarios.

  • 3
    I am not sure that you are right about `pow`. On some implementations, it is pretty well optimized. – Basile Starynkevitch Nov 11 '14 at 08:29
  • 3
    In case of power of two you can simply using shift operations: 1 << N means 2^N. – sfrehse Nov 11 '14 at 08:31
  • 1) Is `b` is an `double` with a whole number value or is `b` an `int`? 2) What type should the result be: `double`, `int`, `unsigned long long, etc? – chux - Reinstate Monica Nov 11 '14 at 08:32
  • @sfrehse `a` is not a power of 2 @legends2k c++11 @chux `b` is `int` and the result is used to perform further calculations with `double`s –  Nov 11 '14 at 08:34
  • 3
    Did you read that it literally performed those function calls? Or just that it had the same effect as doing so? In any case, see http://en.wikipedia.org/wiki/Exponentiation_by_squaring – Mike Seymour Nov 11 '14 at 08:35

2 Answers2

13

There are a number of faster alternatives I've collected over the years that typically rely on a recursive implementation of the function, and bit shifts to handle multiplication when warranted. The following provide functions tailored to integer, float and double. They come with the normal disclaimer: while faster not all possible test have been run and the user should validate input is sane before calling and on return... blah, blah, blah.. But, they are pretty darn useful:

I believe proper attribution goes to Geeks for Geeks Pow(x,n) as pointed out by blue moon. I had long since lost the links.. That looks like them. (minus a tweak or two).

/* Function to calculate x raised to the power y

    Time Complexity: O(n)
    Space Complexity: O(1)
    Algorithmic Paradigm: Divide and conquer.
*/
int power1 (int x, unsigned int y)
{
    if (y == 0)
        return 1;
    else if ((y % 2) == 0)
        return power1 (x, y / 2) * power1 (x, y / 2);
    else
        return x * power1 (x, y / 2) * power1 (x, y / 2);

}

/* Function to calculate x raised to the power y in O(logn)
    Time Complexity of optimized solution: O(logn)
*/
int power2 (int x, unsigned int y)
{
    int temp;
    if (y == 0)
        return 1;

    temp = power2 (x, y / 2);
    if ((y % 2) == 0)
        return temp * temp;
    else
        return x * temp * temp;
}

/* Extended version of power function that can work
for float x and negative y
*/
float powerf (float x, int y)
{
    float temp;
    if (y == 0)
    return 1;
    temp = powerf (x, y / 2);
    if ((y % 2) == 0) {
        return temp * temp;
    } else {
        if (y > 0)
            return x * temp * temp;
        else
            return (temp * temp) / x;
    }
}

/* Extended version of power function that can work
for double x and negative y
*/
double powerd (double x, int y)
{
    double temp;
    if (y == 0)
    return 1;
    temp = powerd (x, y / 2);
    if ((y % 2) == 0) {
        return temp * temp;
    } else {
        if (y > 0)
            return x * temp * temp;
        else
            return (temp * temp) / x;
    }
}
David C. Rankin
  • 81,885
  • 6
  • 58
  • 85
  • Thank you for the detailed answer, however I do not see the bit shifts you were talking about... –  Nov 11 '14 at 08:46
  • 3
    @BlueMoon Thank you for the link. That looks like where they came from originally. On my box they live at `~/dev/src-c/function/math/fn-power.c` and the original links had been lost to the passage of time. I added the attribution. – David C. Rankin Nov 11 '14 at 08:49
  • @Pickle - I must have been thinking of another set where multiplication of less than 10 was broken down into the respective combination of `shifts + adds`. Don't see it in theses. – David C. Rankin Nov 11 '14 at 08:51
  • 1
    `power1()` does not have space complexity `O(1)` but `O(log y)` (note there is the same recursion depth as in the case of `power2()`). And the time complexity of `power1()` is much worse then indicated; it is rather `O(log^2 y)` so it should not be used at all: `power2()` is more optimized version of `power1()` ensuring the same thing is not (on each recursion) calculated twice needlessly making it `O(log y)`. For large `y` this would really matter a lot. (It's true today's compilers may be smart enough to do the optimization for you, but I would not still rely on that.) – mity Mar 05 '18 at 17:22
  • `double` return more accurate results, at least in my case. – vmemmap Aug 11 '22 at 21:25
3

Non-recursive non-floating point answer

Replace uintmax_t/intmax_t with the type of your desire. Overflow not detected.

uintmax_t powjuu(unsigned x, unsigned y) {
  uintmax_t z = 1;
  uintmax_t base = x;
  while (y) {
    if (y & 1) {  // or y%2
      z *= base;
    }
    y >>= 1; // or y /= 2
    base *= base;
  }
  return z;
}

intmax_t powjii(int x, int y) {
  if (y < 0) {
    switch (x) {
      case 0:
        return INTMAX_MAX;
      case 1:
        return 1;
      case -1:
        return y % 2 ? -1 : 1;
    }
    return 0;
  }
  intmax_t z = 1;
  intmax_t base = x;
  while (y) {
    if (y & 1) {
      z *= base;
    }
    y >>= 1;
    base *= base;
  }
  return z;
}
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
  • this returned pow results very fast and I am able to customize it according to my requirement like add a modulus to the pow result – psaraj12 Jun 07 '21 at 06:46
  • @psaraj12 To perform full range `pow(a,b) mod c`, see [Modular exponentiation without range restriction](https://codereview.stackexchange.com/q/187257/29485). – chux - Reinstate Monica Jun 07 '21 at 11:39