Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 3, 4, 5, and 6

Reply
 
Thread Tools Display Modes
  #1  
Old 12-21-2003, 10:59 PM
jrkarthik jrkarthik is offline
Member
 
Join Date: Dec 2003
Location: India
Posts: 2
Unhappy algorithm for calculation of x^(4/3)

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
Reply With Quote
  #2  
Old 12-23-2003, 12:15 PM
TMcCloskey TMcCloskey is offline
Member
 
Join Date: Jun 2002
Location: Stillwater
Posts: 8
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).
Reply With Quote
  #3  
Old 12-23-2003, 10:22 PM
jrkarthik jrkarthik is offline
Member
 
Join Date: Dec 2003
Location: India
Posts: 2
Thumbs up

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..
Reply With Quote
  #4  
Old 01-10-2004, 04:15 AM
jaje jaje is offline
Registered User
 
Join Date: Jul 2002
Posts: 55
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. (^_^)
Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump


All times are GMT -5. The time now is 08:39 PM.


Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2010, Jelsoft Enterprises Ltd.