Recursive power (mathematical operation) function in C

Question

How to make a C recursive exponentiation (power) function better?

I developed a pretty simple (yet decent, imo) power function for the sake of learning. I'm not quite trying to compete with math.h's double pow(double __x, double __y)

double pwr(double _base, double _exp)
{
    if (_exp == 0)
        return 1;
    
    else if (_exp > 0)
        return _base * pwr(_base, _exp - 1);
    else
        return 1 / (_base * pwr(_base, (_exp * -1) - 1));
}

But i would like to know why (if) math.h's would be better. Some community guy also told me it could be further simplified, but didn't tell me how. Is that true?

Note: I'm aware of "divide and conquer" approach, but not sure how to proceed in case _exp is not integer

Answer 1

For large exponents, this code is slower than necessary (it scales linearly in exp, rather than logarithmically as an efficient implementation would). Research Binary Exponentiation algorithm, and Exponentiation by logarithm, which is what your standard library implementation is most likely to use for pow(double, double).

The latter relies on the identity \$\ln{(y^x)} ≡ x \ln{(y)}\$. So \$y^x = \exp(x \ln(y))\$.

The implementation here will never return when exp is not an integer - you'll just alternate between calls either side of zero, until you run out of stack. We need to use the logarithm approach for non-integer powers (and arrange to fail when the base is negative and the power is fractional).

Avoid using reserved identifiers in your own code. This includes names that begin with underscore (_) followed by a lower-case letter. (Not as dangerous as names beginning with _ and an upper-case letter, which are available to the implementation for use as macros, but still not good; as a general rule, never declare your own names beginning with _ and you'll be fine).

Answer 2

For exp not 0, between -1 and 1 it will make a recursive call with a negated exponent.

Furthermore x²ⁿ = (xⁿ)² which can be used to reduce n recursions down to log₂(n) recursions.

For exp = n + q, n int, 0 < q < 1 or exp = n ± q, n int, 0 < q <= 0.5 you can use x^exp = (xⁿ).(x^q). For the integer power n you can do the logarithmic recursion (double, long or BigInteger). For the q fraction you need an approximation sequence. Don't ask me how.

Underscore at the beginning, is extra. Nowadays the majority seems to find this less readable (as not a ubiquitous convention), and disruptive in reading. Less is more here.

Answer 3

pow(double _base, double _exp), like this pwr(), is a challenging function.

Also an accurate implementation is difficult as errors in manipulating _base are heavily amplified by large _exp.

Non-integer power

When _exp is not a whole number value, best to use exp2(_exp * log2(_base)) and avoid log(), exp() to improve accuracy.

Note: _base < 0 will error as, in general, the result of pow(negative_base, some_non_whole_number) is a complex number.

else not needed

double pwr(double _base, double _exp) {
  if (_exp == 0)
    return 1;
  // else if (_exp > 0)
  if (_exp > 0)
    return _base * pwr(_base, _exp - 1);
  // else
  return 1 / (_base * pwr(_base, (_exp * -1) - 1));

Linear recursion depth

Should _exp be a large integer value, code will recurse many times.

Instead, halve the exponent each each recursion for a log recursion depth.

// return _base * pwr(_base, _exp - 1);
double p2 = pwr(_base, _exp/2);
return p2*p2*(_exp/2 ? _base : 1);

This also improves correctness as fewer calculations are done.

---

Somewhat improved code

double pwr_rev2(double base, double expo) {
  // If non-integer expo, use exp2(),log2()
  double expo_ipart;
  double fpart = modf(expo, &expo_ipart);
  if (fpart) {
    return exp2(expo * log2(base));
  }

  if (expo_ipart < 0.0) {
    // Form 1/ result
    return 1.0/pwr_rev2(base, -expo);
  }

  if (expo_ipart == 0.0) {
    return 1.0;
  }

  double p = pwr_rev2(base, expo_ipart / 2.0);
  p *= p;
  if (fmod((expo_ipart, 2)) {
    p *= base;
  }
  return p;
}

Implementation still has short comings with 1.0/pwr_rev2(base, -expo) as the quotient becomes 0.0 when the result should be very small, yet more than zero.

It also recurses too many times when the result p *= p; is infinity. TBD code to fix.