Computing fractional exponents in C

c, exponent, newtons-method

Solution

It seems your task is to calculate

⎛ xN ⎞(aN / aD)
⎜⎼⎼⎼⎼⎟           where xN,xD,aN,aD ∈ ℤ,  xD,aD ≠ 0
⎝ xD ⎠

using only multiplications, divisions, additions, and subtractions, with Newton's method as the suggested method to implement.

The equation we're trying to solve (for y) is

             (aN / aD)
y = (xN / xD)            where y ∈ ℝ

Newton's method finds a root of a function. If we want to use it to solve the above, we substract the right side from the left side, to get a function whose zero gives us the y we want:

                  (aN/aD)
f(y) = y - (xN/xD)        = 0

Not much help. I guess this is as far as you got? The point here is to not form that function just yet, because we don't have a way to calculate a rational power of a rational number!

First, let's decide that aD and xD are both positive. We can do that simply by negating both aN and aD if aD was negative (so sign of aN/aD does not change), and negating both xN and xD if xD was negative. Remember, by definition neither xD or aD is zero. Then, we can simply raise both sides to the aD'th power:

 aD            aN     aN     aN
y   = (xN / xD)   = xN   / xD

We can even eliminate the division by multiplying both sides by the last term:

 aD     aN     aN
y   × xD   = xN

Now, this looks quite promising! The function we get from this is

        aD   aN     aN
f(y) = y   xD   - xN

Newton's method also requires the derivative, which is obviously

f(y)            aD   aN
⎼⎼⎼⎼ = df(y) = y   xD   y / aD
 dy

Newton's method itself relies on iterating

             f(y)
y    = y  - ⎼⎼⎼⎼⎼⎼
 i+1    i    df(y)

If you work out the math, you'll find that the iteration is just

                                 aD
                y[i]      y[i] xN
y[i+1] = y[i] - ⎼⎼⎼⎼ + ⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼⎼
                 aD           aD   aN
                       aD y[i]   xD

You don't need to keep all the y values in memory; it is enough to remember the last one, and stop iterating when their difference is small enough.

You do still have exponentiation above, but now they are integer exponentiation only, i.e.

  aD
xN   = xN × xN × .. × xN
       ╰───────┬───────╯
              aD times

which you can do very simply, for example just by multiplying the argument by itself the desired number of times, e.g. in C,

double ipow(const double base, const int exponent)
{
    double result = 1.0;
    int    i;
    for (i = 0; i < exponent; i++)
        result *= base;
    return result;
}

There are more efficient methods to do integer exponentiation, but the above function should be perfectly acceptable for this.

The final problem is to pick the initial y so that you get convergence. You cannot use 0, because (a power of) y is used as a denominator in the division; you'd get division by zero error. Personally, I'd check whether the result ought to be positive or negative, and smaller than or greater than one in magnitude; two rules overall to pick a safe initial y.

Questions?

Problem

I'm trying to evaluate a^n, where a and n are rational numbers. I don't want to use any predefined functions like `sqrt()` or `pow()` So I'm trying to use Newton's Method to get an approximate solution using this approach: 3^0.2 = 3^(1/5) , so if x = 3^0.2, x^5 = 3. Probably the best way to solve that (without a calculator but still using the basic arithmetic operations) is to use "Newton's method". Newton's method for solving the equation f(x)= 0 is to set up a sequence of numbers xn defined by taking x0 as some initial "guess" and then xn+1= xn- f(xn/f '(xn) where f '(x) is the derivative of f. Posted on physicsforums The problem with that method is that if I want to compute `5.2^0.33333`, I'll need to find the roots for this equation `x^10000 - 5.2^33333 = 0`. I end up with huge numbers, and get `inf` and `nan` errors most of the time. Can someone give me advice on how to solve this problem? Or, can someone provide another algorithm to compute a^n?

Original source