Exponential Function

[chris]

Is there any way I can make this function better or is there any visible problems with it? It seems to act weird with things like 10^10 and above...

Code: Select all

int pow(int base, int exp)
{
   int res = base;
   
   if(exp == 0)
      return 1;
      
   if(exp == 1)
      return base;
   
   for (; exp-1 > 0; exp--) {
      res *= base;
   }      
   
   return res;
   
}

[chris]

Err I guess 10^10 isn't working because I used ints?

[Tim]

Right. Try float, double, or long long.

[sonneveld]

Code: Select all

double pow(double base, int exp)
{
   double res = 1;
   
   while(exp--)
      res *= base;

   return res;
}

- Nick

[chris]

Thanks.

[Schol-R-LEA]

While that will work, it doesn't catch the screw cases (i.e., zero bases, zero exponents, negative exponents), it may have serious round-off errors, and won't be very efficient, either. To really do it right, you need to dig into the double-precision format and manipulate the exponent and mantissa directly. Even without that, though, a few simple checks and modifications can speed things up substantially, as shown in this example (tested for basic function but not for round-off errors):

Code: Select all

#define reciprocal(n) (1 / (n))
#define odd(x) ((x) % 2)
#define even(x) (!odd(x))

double pow(double base, int exp)
{
    double accumulate = 1.0;
    
    if (0.0 == base)    // sanity check - is base zero? 0^exp == 0
    {    
        return 0.0;
    }
    else if (1 == base) // sanity check - is base one? 1^exp == 1 
    {
        return 1.0;
    }
    else if (0 == exp) // sanity check - is exp zero? base^0 == 1
    {
        return 1.0;
    }
    else if (0 > exp)    // base^(-exp)  ==  1 / base^exp
    {
        return reciprocal(pow(base, abs(exp)));
    }
    else 
    
    
    while (0 < exp)
    {
        if (even(exp))
        // by halving even exponents and squaring the base,
        // the number of multiplications needed is reduced
        // turning the algorithm from O(n) time and O(1) space
        // to 0(log n) for time while remaining O(1) for space.
        {
                exp /= 2;
                base *= base;
        }
        else 
        {
        // For odd numbers, multiply the base by base^(exp-1)
        // Since exp - 1 is guaranteed to be even, it will then
        // use the halving formula above to reduce the total 
        // number of multiplications.
                accumulate *= base;
                exp--;
        }
    }
    return accumulate;
}

[chris]

Thanks alot Schol-R-LEA. Comments and everything ;D

[df]

since exp is an integer

replace

Code: Select all

 exp /= 2

with

Code: Select all

exp>>=1

might save you _1_ cycle hehe

[chris]

Code: Select all

        return reciprocal(pow(base, abs(exp)));

One question: Where does the "abs" come from? I didn't notice it before

[Schol-R-LEA]

since exp is an integer

replace

Code: Select all

 exp /= 2

with

Code: Select all

exp>>=1

might save you _1_ cycle hehe

I knew I'd missed something. Actually, it might be more than one cycle, and even if it is only one, well, they do add up...

Actually, the real annoyance to me is that you're doing the division twice, once for the modulo operation in the odd() macro, and then again here. In assembly it wouldn't be a big deal; you'd store the exponent in eax, shift, and jump on the carry flag, but trying to do the equivalent in C is too much work.

In any case, a better approach has occurred to me, at least once you've gotten the exp() and log() functions working:

Code: Select all

double pow(double base, int expt)
{
   return exp((expt * log(base));
}

Writing exp() and log() are left as an exercise (hint: the log() function maps very nicely to the FPU instruction FLDL2E). ;D

[Schol-R-LEA]

Code: Select all

        return reciprocal(pow(base, abs(exp)));

One question: Where does the "abs" come from? I didn't notice it before

Oops. It's exactly the kind of math.h function were supposed to be implementing here, isn't it? OK, then here's a quick and dirty version that should do:

Code: Select all

#define abs(x) (0 < (x)) ? (x) : (-x))

[chris]

Thanks ;D

[Therx]

I've seen it implemented before something like this (off the top of my head)

Code: Select all

double pow(double base, int expt)
{
 if(expt == 0)
  return 1;
 else
  return base * pow(base, expt - 1);
}

With the claim that its multithreaded friendly?

Pete

[Schol-R-LEA]

I've seen it implemented before something like this (off the top of my head)

Code: Select all

double pow(double base, int expt)
{
 if(expt == 0)
  return 1;
 else
  return base * pow(base, expt - 1);
}

With the claim that its multithreaded friendly?

Certainly; this is one of the simplest ways to implement this. There's no reason for it not to be multithread friendly; it's a basic recursive function, with no shared memory state. If you think about it, you'd see that has to be at least serially reentrant, or the recursion wouldn't work properly.

I'd still do it like this, though, to speed it up:

Code: Select all

double pow(double base, int expt)
{
 if(expt == 0)
   return 1;
 else if (expt == 1)
   return base;
 else if (even(base))
  {
   base = pow(base, (expt / 2));
   return base * base;
  }
 else
   return base * pow(base, expt - 1);
}

This is what I based my iterative version on (though I used a slightly different equivalency in the final version); like that one, this is O(n log n) in time, but unlike that version, it is also O(n log n) in space, not O(1). This one is easier to understand, however, though not as simple as yours.