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?
3
- will 'b' always be an integer? if so, just start with 1 and multiply it by a, b timesTom Sirgedas– Tom Sirgedas2010-08-19 05:40:39 +00:00Commented Aug 19, 2010 at 5:40
- a and b are both floating point and will not be natural numbersymihere– ymihere2010-08-19 05:46:12 +00:00Commented Aug 19, 2010 at 5:46
- do you have sqrt() available?Tom Sirgedas– Tom Sirgedas2010-08-19 06:01:06 +00:00Commented Aug 19, 2010 at 6:01
Add a comment |
3 Answers
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 
5 Comments
ymihere
thanks alot. this is precisely was what i was looking for. Out of interest: can you give me any background to that algorithm?
Tom Sirgedas
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.
Drag-On
Note that, if your sqrt function has the same restrictions as std::sqrt, then this won't work for negative bases.
Sean Tashlik
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)).
Eduard Grigoryev
Btw, you can implement sqrt as in stackoverflow.com/a/3047531/8990391.
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.
4 Comments
ymihere
i have a working ln function. However, for the Taylor series I need powers again.
Jim Lewis
@ymihere: The Taylor series expansion only contains integer exponents, which can be reduced to multiplication.
Tom Sirgedas
@ymihere: do you have exp() available? if so, this solution is best!
ymihere
@tom: i don't have exp. actually, that's what I am working around.
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; } } 