The most efficient way of implementing pow() function in floating point

Question

I am trying to implement my own version of pow() and sqrt() function as my custom library doesn't have pow()/sqrt() floating point support.

Can anyone help?

Answer 1

Sure - it's easy if you have exponential and natural log functions.

Since y = x^n, you can take the natural log of both sides:

ln(y) = n*ln(x)

Then taking the exponential of both sides gives you what you want:

y = exp(n*ln(x))

If you want something better, the best place I know to look is Abramowitz and Stegun.

Answer 2

Yes, Sun can (Oracle now, I guess):

fdlibm, the "freely distributable math library", has sqrt and pow, along with many other math functions.

They're fairly high-tech implementations, though, and of course nothing is ever the "most efficient" implementation of something like this. Are you after source code to get it done, or are you really not so much looking for pow and sqrt, but actually looking for an education in floating-point algorithms programming?

Answer 3

Note that if your instruction set has an instruction for square root or power, you'll be much better off using that. The x87 floating point instructions, for example, have an instruction fsqrt, and the SSE2 additions include another instruction sqrtsd, which are probably going to be much faster than most solutions written in C. In fact, atleast gcc uses the two instructions when compilation takes place on an x86 machine.

For power, however, things get somewhat murky. There's an instruction in the x87 floating point instruction set that can be used to calculate n*log2(n), namely fyl2x. Another instruction, fldl2e, stores log2(e) in the floating point stack. You might want to give these a look.

You might also want to take a look at how individual C libraries do this. dietlibc, for example, simply uses fsqrt:

sqrt:
    fldl 4(%esp)
    fsqrt
    ret

glibc uses Sun's implementation for machines where a hardware square root instruction is not available (under sysdeps/ieee754/flt-32/e-sqrtf.c), and uses fsqrt on the x86 instruction set (though gcc can be instructed to instead use the sqrtsd instruction.)

Answer 4

Square root is properly implemented with an iterative Newtons method.

Answer 5

double ipow(int base, int exp)
{

bool flag=0;
if(exp<0) {flag=1;exp*=-1;}
int result = 1;
while (exp)
{
    if (exp & 1)
        result *= base;
    exp >>= 1;
    base *= base;
}
if(flag==0)
return result;
else
return (1.0/result);
}
//most suitable way to implement power function for integer to power integer

Answer 6

For calculating the square root of a float in C I'd recommend using fsqrt if you target x86. You can use such ASM instruction with:

asm("fsqrt" : "+t"(myfloat));

For GCC or

asm {

fstp myfloat

fsqrt

fldp myfloat

}

Or something like that for Visual Studio.

For implementing pow, using a big switch statement like the one at upitasoft.com/link/powLUT.h should do. It can cause some cache problems but if you keep it like that it shouldn't be a problem, just limit the range (note, you can still optimize the code I provided).

If you want to support floating point powers, is way harder... You can try using the natural logarithm and exponential functions, such as:

float result = exp(number * log(power));

But usually it is slow and/or imprecise.

Hope I helped.

Answer 7

The fastest way I can think of doing a pow() would be along these lines (note, this is pretty complicated):

//raise x^y
double pow(double x, int y) {
    int power;
    map<int, double> powers;
    for (power = 1; power < y; power *= 2, x *= x)
        powers.insert(power, x);

    while (power > y) {
        //figure out how to get there
        map<int, double>::iterator p = powers.lower_bound(power - y);
        //p is an iterator that points to the biggest power we have that doesn't go over power - y
        power -= p->first;
        x /= p->second;
    }

    return x;
}

I have no idea about how to implement a decimal power. My best guess would be to use logarithms.

Edit: I'm attempting a logarithmic solution (based on y), as opposed to a linear solution, which you propose. Let me work this out and edit it, because I know it works.

Edit 2: Hehe, my bad. power *= 2 instead of power++