Computing fractional exponents in C

Question

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?

Answer 1

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?

Answer 2

You can use the generalized binomial theorem. Substitute y=1 and x=a-1. You would want to truncate the infinite series after enough terms, based on the desired accuracy. To be able to link number of terms to accuracy, you would need to ensure that the x^r terms are decreasing in absolute value. So, depending on the value of a and n, you should apply the formula to compute one of a^n and a^(-n) and use that to get your desired result.

Answer 3

A solution for raising an integer number to a power is:

int poweri (int x, unsigned int y)
{
    int temp;
    if (y == 0)
        return 1;

    temp = poweri (x, y / 2);
    if ((y % 2) == 0)
        return temp * temp;
    else
        return x * temp * temp;
}

However, the square root doesn't provide as clean of a closed solution. There is a good bit of background to be found at wikipedia-square root and at Wolfram Mathworks Square Root Algorithms Both provide several methods that will meet your needs, you just have to choose the one that fits your purpose.

With slight modification, this routine from wikipedia (modified to return the square root and refine accuracy) returns a surprisingly accurate square root. Yes, there will be howls about the use of a union, and it is only valid where integer and float storage are equivalent, but if you are hacking your own square root, this is relatively efficient:

float sqrt_f (float x)
{
        float xhalf = 0.5f*x;
        union
        {
            float x;
            int i;
        } u;
        u.x = x;
        u.i = 0x5f3759df - (u.i >> 1);
        /* The next line can be repeated any number of times to increase accuracy */
        // u.x = u.x * (1.5f - xhalf * u.x * u.x);
        int i = 10;
        while (i--)
            u.x *= 1.5f - xhalf * u.x * u.x;

        return 1.0f / u.x;
}