Floating point exponentiation without power-function

Question

Currently I have to work in an environment where the power-operator is bugged. Can anyone think of a method temporarily work around this bug and compute a^b (both floating point) without a power function or operator?

Answer 1

You can use the identity a^b = e^{(b log a)}, then all the calculations are relative to the same base e = 2.71828...

Now you have to implement f(x) = ln(x), and g(x) = e^x. The fast, low precision method would be to use lookup tables for f(x) and g(x). Maybe that's good enough for your purposes. If not, you can use the Taylor series expansions to express ln(x) and e^x in terms of multiplication and addition.

Answer 2

if you have sqrt() available:

double sqr( double x ) { return x * x; }
// meaning of 'precision': the returned answer should be base^x, where
//                         x is in [power-precision/2,power+precision/2]
double mypow( double base, double power, double precision )
{   
   if ( power < 0 ) return 1 / mypow( base, -power, precision );
   if ( power >= 10 ) return sqr( mypow( base, power/2, precision/2 ) );
   if ( power >= 1 ) return base * mypow( base, power-1, precision );
   if ( precision >= 1 ) return sqrt( base );
   return sqrt( mypow( base, power*2, precision*2 ) );
}
double mypow( double base, double power ) { return mypow( base, power, .000001 ); }

test code:

void main()
{
   cout.precision( 12 );
   cout << mypow( 2.7, 1.23456 ) << endl;
   cout << pow  ( 2.7, 1.23456 ) << endl;
   cout << mypow( 1.001, 1000.7 ) << endl;
   cout << pow  ( 1.001, 1000.7 ) << endl;
   cout << mypow( .3, -10.7 ) << endl;
   cout << pow  ( .3, -10.7 ) << endl;
   cout << mypow( 100000, .00001 ) << endl;
   cout << pow  ( 100000, .00001 ) << endl;
   cout << mypow( 100000, .0000001 ) << endl;
   cout << pow  ( 100000, .0000001 ) << endl;
}

outputs:

3.40835049344
3.40835206431
2.71882549461
2.71882549383
393371.348073
393371.212573
1.00011529225
1.00011513588
1.00000548981
1.00000115129

Answer 3

given that you can use sqrt, this simple recursive algorithm works:

Suppose that we're calculating aˆb. The way the algorithm works is by doing Fast Exponentiation on the exponent until we hit the fractional part, once in the fractional part, do a modified binary search, until we're close enough to the fractional part.

double EPS = 0.0001;

double exponentiation(double base, double exp){
  if(exp >= 1){
    double temp = exponentiation(base, exp / 2);
    return temp * temp;
  } else{
    double low = 0;
    double high = 1.0;

    double sqr = sqrt(base);
    double acc = sqr;    
    double mid = high / 2;

    while(abs(mid - exp) > EPS){
      sqr = sqrt(sqr);

      if (mid <= exp) {
          low = mid;
          acc *= sqr;
      } else{
          high = mid;
          acc *= (1/sqr);
      }

      mid = (low + high) / 2;
    }

    return acc;
  }
}