10

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?

ymihere
  • 381
  • 1
  • 4
  • 11

3 Answers3

26

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
Tom Sirgedas
  • 3,208
  • 20
  • 17
  • thanks alot. this is precisely was what i was looking for. Out of interest: can you give me any background to that algorithm? – ymihere Aug 19 '10 at 07:40
  • 7
    Sure, the basic idea is that x^.5 = sqrt(x), x^.25 = sqrt(sqrt(x)), x^.125 = sqrt(sqrt(sqrt(x))), etc. With these building blocks, we can say x^.625 = (x^.5)*(x^.125). We can't express, say, x^.3 exactly, but we can get arbitrarily close. I implemented this a little differently, but it uses the same concept. – Tom Sirgedas Aug 19 '10 at 15:28
  • Note that, if your sqrt function has the same restrictions as std::sqrt, then this won't work for negative bases. – Drag-On Aug 11 '19 at 14:59
  • I like your idea, but I guess it can be done in a better complexity in both time and space. The iterative version can save space up to O(1) and time can be improved up to O(log(n)). – Sean Tashlik Sep 13 '20 at 12:58
  • Btw, you can implement sqrt as in https://stackoverflow.com/a/3047531/8990391. – Eduard Grigoryev Mar 28 '21 at 12:57
10

You can use the identity ab = 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.

Jim Lewis
  • 43,505
  • 7
  • 82
  • 96
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;
  }
}
eapolinario
  • 105
  • 7