Newton Raphson with SSE2 - can someone explain me these 3 lines

Question

I'm reading this document: http://software.intel.com/en-us/articles/interactive-ray-tracing

and I stumbled upon these three lines of code:

The SIMD version is already quite a bit faster, but we can do better. Intel has added a fast 1/sqrt(x) function to the SSE2 instruction set. The only drawback is that its precision is limited. We need the precision, so we refine it using Newton-Rhapson:

 __m128 nr = _mm_rsqrt_ps( x ); 
 __m128 muls = _mm_mul_ps( _mm_mul_ps( x, nr ), nr ); 
 result = _mm_mul_ps( _mm_mul_ps( half, nr ), _mm_sub_ps( three, muls ) ); 

This code assumes the existence of a __m128 variable named 'half' (four times 0.5f) and a variable 'three' (four times 3.0f).

I know how to use Newton Raphson to calculate a function's zero and I know how to use it to calculate the square root of a number but I just can't see how this code performs it.

Can someone explain it to me please?

Splanchnic

Answer 1

Given the Newton iteration , it should be quite straight forward to see this in the source code.

 __m128 nr   = _mm_rsqrt_ps( x );                  // The initial approximation y_0
 __m128 muls = _mm_mul_ps( _mm_mul_ps( x, nr ), nr ); // muls = x*nr*nr == x(y_n)^2
 result = _mm_mul_ps(
               _mm_sub_ps( three, muls )    // this is 3.0 - mul;
   /*multiplied by */ __mm_mul_ps(half,nr)  // y_0 / 2 or y_0 * 0.5
 );

And to be precise, this algorithm is for the inverse square root.

Note that this still doesn't give fully a fully accurate result. rsqrtps with a NR iteration gives almost 23 bits of accuracy, vs. sqrtps's 24 bits with correct rounding for the last bit.

The limited accuracy is an issue if you want to truncate the result to integer. (int)4.99999 is 4. Also, watch out for the x == 0.0 case if using sqrt(x) ~= x * sqrt(x), because 0 * +Inf = NaN.

Hypoacidity

When truncating to integer, do you think it would be viable as a final step to add a value which has the same exponent as the result but only the lowest one (or two?) bits set in the significand? This is of course under the condition that the least significant figure is always lower than the one's position. – Harsh

It's application dependent. The point is that when using an iterative approach sqrt(n*n) == n doesn't always hold. This can't be "fixed" arbitrarily -- as sqrt(n*n - epsilon) == n may lead to disaster. – Hypoacidity

Answer 2

To compute the inverse square root of a, Newton's method is applied to the equation 0=f(x)=a-x^(-2) with derivative f'(x)=2*x^(-3) and thus the iteration step

N(x) = x - f(x)/f'(x) = x - (a*x^3-x)/2 
     = x/2 * (3 - a*x^2)

This division-free method has -- in contrast to the globally converging Heron's method -- a limited region of convergence, so you need an already good approximation of the inverse square root to get a better approximation.

Kantianism