7

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?

1
  • I expect the answer depends strongly on which Instruction Set you're targeting. It's hard to answer for the general case, except at a very broad level. 2 days ago

7 Answers 7

12

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?

7
  • @Steve:- Firstly to find out a better source code with good precision and secondly to actually educate myself on this topic.
    – Viks
    Apr 19, 2010 at 0:18
  • 1
    For the latter I think you really need a good book - there's only so much you can achieve staring at source code, especially when it's mostly made of magic numbers and corner-cases. I'm afraid it's not an area I'm strong on, though: I recommend fdlibm solely because I know that it works tolerably well, not because I actually understand the source. Apr 19, 2010 at 0:27
  • 1
    Most of us can do something halfway clever with a Taylor series to implement a function, but finding a really fast way to compute such a function is going to take specialized knowledge, and may vary depending on platform. You're probably better off just taking somebody else's source and running. Apr 19, 2010 at 20:56
  • 2
    fdlibm is an excellent starting point. It's not the fastest implementation on any given platform, but it is portable and good enough for most purposes. If you want the most efficient implementation on a given platform, you will need to spend a few years learning about the micro-architecture of your target platform, then a few more years learning low-level numerics -- or hire someone who has already done so =) Apr 19, 2010 at 21:02
  • 1
    The code for the approximation with a taylor series expansions are way more slower than other approximations...
    – Quonux
    Apr 4, 2013 at 9:21
12

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.

8
  • 2
    And if you don't have exponential and log functions, you can use a Taylor series approximation. But I'm not sure if that is the "most efficient" either...
    – rlbond
    Apr 18, 2010 at 23:54
  • For x > 0, anyway. With negative input not so much. Apr 18, 2010 at 23:55
  • Even for x > 0, you'll lose precision doing it this way, compared to a sophisticated specialised implementation such as that found in fdlibm. Apr 19, 2010 at 0:09
  • 2
    rlbond: a Taylor series is a bad idea. Chebyshev polynomials are probably better.
    – Gabe
    Apr 19, 2010 at 1:22
  • 5
    More to the point, exp(y log x) delivers much lower accuracy than a does a proper implementation of pow(x,y). That may be acceptable for your purposes, or it may not. Apr 19, 2010 at 21:00
5

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.)

2

Square root is properly implemented with an iterative Newtons method.

1
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
0

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.

-2

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++

3
  • The practicality of this algorithm aside, as written it doesn't work. Consider a simple case, x = 2, y = 3. power = 1, x = 4. power = 2, x = 8. power = 3, x = 16. power == y, so the first loop stops. power is still == y, so the second loop doesn't do anything. x == 16, but 2^3 == 8. Further, the only time the second loop would ever do anything is if y < 0; I'll leave why that's a disaster as an exercise for the reader. Apr 19, 2010 at 1:12
  • Erp, my summary was wrong. You're squaring x each iteration, rather than multiplying by your original value. x -> xx -> xxxx when you want x -> xx -> xx*x. But like I said, your second loop [and by extension, the entire map you build] is completely pointless because your first loop will end when power == y unless y is negative. Apr 19, 2010 at 1:36
  • 2
    construct/populate/use/destruct a map in a frequently called function is a questionable activity.
    – EvilTeach
    Jul 12, 2013 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.