Self made pow () C ++

Question

I read through How can I write a power function myself? , and the answer given by dan04 attracted my attention mainly because I am not sure about the answer given by fortran, but I accepted it and executed it:

#include <iostream> using namespace std; float pow(float base, float ex){ // power of 0 if (ex == 0){ return 1; // negative exponenet }else if( ex < 0){ return 1 / pow(base, -ex); // even exponenet }else if ((int)ex % 2 == 0){ float half_pow = pow(base, ex/2); return half_pow * half_pow; //integer exponenet }else{ return base * pow(base, ex - 1); } } int main(){ for (int ii = 0; ii< 10; ii++){\ cout << "pow(" << ii << ".5) = " << pow(ii, .5) << endl; cout << "pow(" << ii << ",2) = " << pow(ii, 2) << endl; cout << "pow(" << ii << ",3) = " << pow(ii, 3) << endl; } } 

although I’m not sure that I translated this right because all calls giving .5 as an exponent return 0. The answer says he might need log2 (x) based on a^b = 2^(b * log2(a)) , but I'm unsure of putting it as I am not sure where to express it, or if I even think about this right.

NOTE. I know that this can be defined in the math library, but I don’t need all the additional costs of the entire math library for several functions.

EDIT: Does anyone know a floating point implementation for fractional metrics? (I saw a double implementation, but it used a register trick, and I need a floating point, and adding a library just to do the trick would be better for me, including just the math library)

Answer 1

Below are links to real powf implementations.

I expect that simpler solutions will not have enough accuracy as a result or the InF and NaN parameters will not be processed.

• http://opensource.apple.com/source/Libm/Libm-2026/Source/Intel/expf_logf_powf.c
• http://opensource.apple.com/source/Libm/Libm-315/Source/ARM/powf.c

Answer 2

I reviewed this article here , which describes how to approximate an exponential function for double precision. After a little research on Wikipedia on the single-precision floating-point representation, I developed equivalent algorithms. They only implemented the exp function, so I found the inverse function for the log, and then just did

  POW(a, b) = EXP(LOG(a) * b). 

compiling this gcc4.6.2 gives the pow function almost 4 times faster than the standard library implementation (compiling with O2).

Note: the code for EXP is copied almost verbatim from the paper I read and the LOG function is copied from here .

Here is the relevant code:

  #define EXP_A 184 #define EXP_C 16249 float EXP(float y) { union { float d; struct { #ifdef LITTLE_ENDIAN short j, i; #else short i, j; #endif } n; } eco; eco.ni = EXP_A*(y) + (EXP_C); eco.nj = 0; return eco.d; } float LOG(float y) { int * nTemp = (int*)&y; y = (*nTemp) >> 16; return (y - EXP_C) / EXP_A; } float POW(float b, float p) { return EXP(LOG(b) * p); } 

There is still some kind of optimization that you can do here, or maybe it's good enough. This is a rough approximation, but if you were satisfied with the errors introduced using the double representation, I assume that it would be satisfactory.

Answer 3

I think you could try to solve it using the Taylor series. Check this out. http://en.wikipedia.org/wiki/Taylor_series

With the Taylor series, you can solve any calculation difficult to solve, such as 3 ^ 3.8, using already known results, such as 3 ^ 4. In this case, you have 3 ^ 4 = 81, therefore

3 ^ 3.8 = 81 + 3.8 * 3 (3.8 - 4) + .. + .. etc. depend on how big your n is, you get a closer solution to your problem.

Answer 4

I think the algorithm you are looking for may be nth root '. Under the initial assumption 1 (for k == 0):

 #include <iostream> using namespace std; float pow(float base, float ex); float nth_root(float A, int n) { const int K = 6; float x[K] = {1}; for (int k = 0; k < K - 1; k++) x[k + 1] = (1.0 / n) * ((n - 1) * x[k] + A / pow(x[k], n - 1)); return x[K-1]; } float pow(float base, float ex){ if (base == 0) return 0; // power of 0 if (ex == 0){ return 1; // negative exponenet }else if( ex < 0){ return 1 / pow(base, -ex); // fractional exponent }else if (ex > 0 && ex < 1){ return nth_root(base, 1/ex); }else if ((int)ex % 2 == 0){ float half_pow = pow(base, ex/2); return half_pow * half_pow; //integer exponenet }else{ return base * pow(base, ex - 1); } } int main_pow(int, char **){ for (int ii = 0; ii< 10; ii++){\ cout << "pow(" << ii << ", .5) = " << pow(ii, .5) << endl; cout << "pow(" << ii << ", 2) = " << pow(ii, 2) << endl; cout << "pow(" << ii << ", 3) = " << pow(ii, 3) << endl; } return 0; } 

Test:

 pow(0, .5) = 0.03125 pow(0, 2) = 0 pow(0, 3) = 0 pow(1, .5) = 1 pow(1, 2) = 1 pow(1, 3) = 1 pow(2, .5) = 1.41421 pow(2, 2) = 4 pow(2, 3) = 8 pow(3, .5) = 1.73205 pow(3, 2) = 9 pow(3, 3) = 27 pow(4, .5) = 2 pow(4, 2) = 16 pow(4, 3) = 64 pow(5, .5) = 2.23607 pow(5, 2) = 25 pow(5, 3) = 125 pow(6, .5) = 2.44949 pow(6, 2) = 36 pow(6, 3) = 216 pow(7, .5) = 2.64575 pow(7, 2) = 49 pow(7, 3) = 343 pow(8, .5) = 2.82843 pow(8, 2) = 64 pow(8, 3) = 512 pow(9, .5) = 3 pow(9, 2) = 81 pow(9, 3) = 729 

Answer 5

My friend and I faced a similar problem while we are working on the OpenGL project, and math.h is not enough in some cases. Our instructor also had the same problem, and he told us to divide the power into whole and floating parts. For example, if you want to calculate x ^ 11.5, you can calculate sqrt (x ^ 115, 10), which can lead to a more accurate result.