Inverse sqrt for fixed point

Question

I am looking for the best inverse square root algorithm for fixed point 16.16 numbers. The code below is what I have so far(but basically it takes the square root and divides by the original number, and I would like to get the inverse square root without a division). If it changes anything, the code will be compiled for armv5te.

uint32_t INVSQRT(uint32_t n)
{
    uint64_t op, res, one;
    op = ((uint64_t)n<<16);
    res = 0;
    one = (uint64_t)1 << 46;
    while (one > op) one >>= 2;
    while (one != 0)
    {
        if (op >= res + one)
        {
            op -= (res + one);
            res +=  (one<<1);
        }
        res >>= 1;
        one >>= 2;
    }
    res<<=16;
    res /= n;
    return(res);
}

Answer 1

The trick is to apply Newton's method to the problem x - 1/y^2 = 0. So, given x, solve for y using an iterative scheme.

Y_(n+1) = y_n * (3 - x*y_n^2)/2

The divide by 2 is just a bit shift, or at worst, a multiply by 0.5. This scheme converges to y=1/sqrt(x), exactly as requested, and without any true divides at all.

The only problem is that you need a decent starting value for y. As I recall there are limits on the estimate y for the iterations to converge.