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.