views:

633

answers:

2

I have a application which needs to raise a number to a fractional power. The target platform is an FPGA and I can get estimates on an FPU size for it, but I need an algorithm for raising a number to a fractional power just for a feasibility study. I'm assuming floating point as a worst case, I expect in practice we will be able to use short cuts, but for now I want to show that worst case can be implemented on our part.

Thought I'd ask here and see if there were any common methods that I could checkout. I know there are software methods of doing this, I want just a reasonably efficient algorithm to start with. I'll worry about the FPGA implementation.

+3  A: 

Is your range of inputs arbitrary, or known within a certain range?

in either case, xm = exp(m log x), so if you can create functions to evaluate exp(x) and log(x) and have a multiply, you're probably all set.

You'll have to figure out how you want to handle nonpositive values of x.

(hint for log(x): if this is IEEE-754 floating point, shift the significand if necessary until you end up with a range of numbers between 2k and 2k+1 for some value K. This lets you deal with a 2:1 range which isn't too hard to approximate with polynomials. Then you just have a small number of possibilities to handle for the exponent and the shift number.

Corresponding hint for exp(x): write x = k+b where 0 <= b < 1 and k is an integer. Then exp(x) = exp(k)*exp(b); b has limited range and k has a limited number of discrete possibilities.)

(hint #2: the numbers probably work out better for xm = g(m f(x)) where f(x) = log2x and g(x) = 2x.)

Jason S
I'm assuming arbitrary for now, I have to read some more of the system's theory behind the calculation were performing. Right now were just trying to propose a solution to the customer. I just want to analytically prove it CAN be done on our system.
NoMoreZealots
+1. I worked on a weird runtime a while back that lacked a full FP library but provided frexp/ldexp routines to separate and recombine the exponent and mantissa; from this, building the standard FP functions was fairly straightforward.
Jeffrey Hantin
Actually for one case in particular it looks like the exponent is going to be > 1. I'm not sure that holds true for all cases.
NoMoreZealots
+2  A: 

As Jason S said, this is done using the identity xm = exp(m log x). In practice though, you will have to deal with truncation errors. I think that the way this is usually done is

  1. Use this identity: xm = (2n * x / 2n)m = 2nm * (x/2n)m and find an integer n such that 1 <= x/2n < 2.
  2. Calculate t = log2(x / 2n). This can be done either using a sufficiently high degree Taylor expansion, or with good ol' Newton-Raphson. You have to make sure that the maximum error over the interval [1, 2[ is not too large for you.
  3. Calculate u = nm + tm.
  4. Our goal is now to calculate 2u. Use the fact that 2u = 2v * 2u-v and find an integer v such that 0 <= u-v < 1.
  5. Calculate w = 2u-v, again using our friend Taylor or Newton-Raphson. The interval of interest is 0 <= u-v < 1.
  6. Your answer is now 2v * w.
erikkallen
How exactly do you use newton raphson to calc a log function? I can find divison algo using it, but there is too much "noise" on google to find how to calc log using it.
NoMoreZealots
I think I remember from my calculation theory that it's possible, but I might be wrong. I'm not a mathematician.
erikkallen