The most efficient way to implement pow () in a floating point

Question

I am trying to implement my own version of the pow () and sqrt () functions, as my user library does not support pow () / sqrt () floating point support.

Can anyone help?

Answer 1

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

fdlibm , the "free math library", sqrt and pow , along with many other math functions.

These are pretty high-tech implementations, and of course, nothing is ever the “most effective” implementation of something like that. Are you following the source code to do this, or are you really not so much looking for pow and sqrt , but are you really looking for education in floating point programming?

Answer 2

Of course, this is easy if you have the exponential and natural functions of a journal.

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

 ln(y) = n*ln(x) 

Then, taking the exponent of both sides, you get what you want:

 y = exp(n*ln(x)) 

If you need something better, the best place I know is Abramowitz and Stegun .

Answer 3

Please note that if your instruction set has instructions for square root or power, you will be much better off using this. For example, x71 floating point instructions have the fsqrt command, and SSE2 add-ons include another sqrtsd command, which is likely to be much faster than most solutions written in C. In fact, at least gcc uses two when compiling is done on x86 machine.

For the authorities, however, everything is becoming somewhat muddy. There is an instruction in the x87 floating-point instruction set, which can be used to calculate n * log2 (n), namely fyl2x . Another fldl2e command stores fldl2e (e) on a floating point stack. You can look at them.

You can also see how separate C. dietlibc libraries do this, for example, just use fsqrt :

 sqrt: fldl 4(%esp) fsqrt ret 

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

Answer 4

The square root is correctly implemented using the Newtons iterative 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

To calculate the square root of a float in C, I would recommend using fsqrt if you are targeting x86. You can use such an ASM statement with:

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

For gcc or

 asm { fstp myfloat fsqrt fldp myfloat 

}

Or something similar for Visual Studio.

To implement pow, use the big switch command, such as upitasoft.com/link/powLUT.h . This may cause some cache problems, but if you keep it that way it shouldn't be a problem, just limit the range (note, you can still optimize the code that I provided).

If you want to support floating point capabilities, this is much more complicated ... You can try using the natural logarithm and exponential functions, such as:

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

But usually it’s slow and / or inaccurate.

Hope I helped.

Answer 7

The fastest way that I can think about what pow () will do is 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 do not know how to realize decimal power. It would be best to use logarithms.

Edit: I am trying to use a logarithmic solution (based on y), not the linear solution that you are proposing. Let me work through this and edit, because I know that it works.

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