self made pow() c++

Question

I was reading through How can I write a power function myself? and the answer given by dan04 caught my attention mainly because I am not sure about the answer given by fortran, but I took that and implemented this:

#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;
    }
}

though I am not sure if I translated this right because all of the calls giving .5 as the exponent return 0. In the answer it states that it might need a log2(x) based on a^b = 2^(b * log2(a)), but I am unsure about putting that in as I am unsure where to put it, or if I am even thinking about this right.

NOTE: I know that this might be defined in a math library, but I don't need all the added expense of an entire math library for a few functions.

EDIT: does anyone know a floating-point implementation for fractional exponents? (I have seen a double implementation, but that was using a trick with registers, and I need floating-point, and adding a library just to do a trick I would be better off just including the math library)

Answer 1

I have looked at this paper here which describes how to approximate the exponential function for double precision. After a little research on Wikipedia about single precision floating point representation I have worked out the equivalent algorithms. They only implemented the exp function, so I found an inverse function for the log and then simply did

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

compiling this gcc4.6.2 yields a pow function almost 4 times faster than the standard library's 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.n.i = EXP_A*(y) + (EXP_C);
      eco.n.j = 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 optimization you can do here, or perhaps that is good enough. This is a rough approximation but if you would have been satisfied with the errors introduced using the double representation, I imagine this will be satisfactory.

Answer 2

I think the algorithm you're looking for could be 'nth root'. With an initial guess of 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 3

I think that you could try to solve it by using the Taylor's series, check this. http://en.wikipedia.org/wiki/Taylor_series

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

3^3.8 = 81 + 3.8*3( 3.8 - 4) +..+.. and so on depend on how big is your n you will get the closer solution of your problem.

Answer 4

I and my friend faced similar problem while we're on an OpenGL project and math.h didn't suffice in some cases. Our instructor also had the same problem and he told us to seperate power to integer and floating parts. For example, if you are to calculate x^11.5 you may calculate sqrt(x^115, 10) which may result more accurate result.