jrkarthik
12-21-2003, 10:59 PM
Hi,
I'm seeking algorithms for computation of x^(4/3). This is for implementation in assembly language. An effective way of implementing Newton Raphson technique or any other algorithm is what i seek. Can anybody help..?
jrk
TMcCloskey
12-23-2003, 12:15 PM
I would recommend approaching this problem as finding a solution to y= x * x^(1/3), ie, let Y = X * CBRT(X) where CBRT is the cube root function.
I say this as you should be able to find proposed solutions to CBRT in the public domain of various approaches/algorithms.
One approach would be:
Let
X = 2^(3p-q)*u where q = 0, 1, or 2 and 1/2 <= u <1.
Then
CBRT(X) = 2^p * CubeRoot(1/2)^q * CubeRoot(u)
This can be summarize as:
(a) Range Reduction: reduce range of X by decomposing X into mantissa (1/2 <= u < 1) and exponent (e = 3p-q). This is relatively easy to do with IEEE real numbers, especially in assembly.
(b) Initial estimate:
- First, obtain an approximation to CubeRoot(u). Call our approximation function "g(u)". Since the range of u is now restricted to 1/2<=u<1, we can precalculate an approximation for g = CubeRoot(u) to a desired precision. This precision does not need to be high since we will iterate to precision in following steps. The number of these iterations, however, will determine the precision required for this initial approximation. This approximation can be via polynomial, rational, or other suitable approximation where we know max error behaviour.
- Second, compute exponent remainders. If we compute p = (e+2)\3 (where "\" represent integer division), then let q = 3p-q. We now have, via simple integer math, the required p and, more importantly, a value of q that is 0, 1, or 2. Thus, there are only two states of q that are none zero. We calculate and store the two values of:
q=1: C(1) = CubeRoot(1/2)^1 = .793700525984....
q=2: C(2) = CubeRoot(1/2)^2 = .629960524947....
- third, compute the initial estimate :
Yo = 2^p * c *g(u) ' use c = 1 if q = 0, c = C(q) otherwise
(c) Now that you have inital estimate to a known precision,you can compute the final answer to desired precision with a finite (ie, predetermined) number of iterations.
For the iteration scheme, consider the follwing
We want to solve the nonlinear equation
F(y) = y^3 - x = 0
Utilize Newton Method
y(i+1) = y(i) - F(y(i))/F'(y(i))
y(i+1) = (1/3) *(2y(i) + x/y(i)^2)
Note, in our iteration scheme, "(1/3)" would be represented as a constant, the "y^2" would be "y*y", the "2y" would be "y+y", etc, to save clock cycles wherever possible (for speed).
Finally, multiply your final answer by "x" to obtain desired function x^(4/3).
jrkarthik
12-23-2003, 10:22 PM
Thanks a lot TMcCloskey! By the way,the implementation is on a DSP.. so no worries abt multiplication! And splitting x^(4/3) into x*cbrt(X) is what i'm doing too. The initialisation is done by polynomial approximation method and i'm settling for the accuracy given by 3 iterations(my inital estimate is good enough..).
I'll try the scale reduction technique and check if its faster..but the input data is audio data and so the range of expected values is known.A look-up table has been used as a considerable part of the data are in reasonable range to have a look-up table for them.;) So..the chances of this technique being faster(For this data set-audio) are pretty low. Anyway i'll try...
thanks again..
If you wish, you can also use the following iteration for computing the cube root (using the same initial estimates mentioned in TMcCloskey's post):
y(i+1) = (2*x*y(i)+y(i)^4)/(x+2*y(i)^3)
This converges in less steps than Newton's method.
Of course, you should do a timing comparison to determine which algorithm is more appropriate for your situation.
Jan M. (^_^)