Is there an fast algorithm computing powers multiple a one-half?

Question

I'm writing a program for solving of planar restricted three-body problem. Its equations are below. This function computes derivatives of position and velocity and write them to array.

valarray<double> force(double t, valarray<double> r)
{
    valarray<double> f(dim);
    valarray<double>r0(r-rb0);
    valarray<double>r1(r-rb1);      

    f[0]=   2 * r[1] + r[2] - (1 - mu)*r0[2]/norm3(r0) - mu*r1[2]/norm3(r1);
    f[1]= - 2 * r[0] + r[3] - mu*r0[3]/norm3(r0) - mu*r1[3]/norm3(r1);
    f[2] = r[0];
    f[3] = r[1];
    return f;
}

double norm3(valarray<double> x)
{
    return pow(x[2]*x[2]+x[3]*x[3],1.5);
}

So I have to compute square of position vector and then raise it to power of 3/2. I think that these operations take a big part of computation time.

Now I use pow function of math.h. Is there another faster algorithm for computation this power? I tried to use fast inverse square root (and cube it later) but it gives too imprecise value for my purposes and works longer (perhaps because of cubing).

Thanks!

Answer 1

A simple approach might be try x*sqrt(x), but benchmark it to be sure.

double norm3(valarray<double> x)
{
    double result=x[2]*x[2]+x[3]*x[3];
    result=result * sqrt(result);
    return result;
}