Newton Raphson与SSE2 - 有人可以解释我这3行

Question

我正在阅读这份文件:http://software.intel.com/en-us/articles/interactive-ray-tracing

我偶然发现了这三行代码:

SIMD版本已经快了很多,但我们可以做得更好.英特尔为SSE2指令集添加了快速1/sqrt(x)函数.唯一的缺点是它的精度有限.我们需要精度,所以我们使用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 ) ); 

此代码假定存在名为"half"(四次0.5f)和变量"three"(四次3.0f)的__m128变量.

我知道如何使用牛顿拉夫森计算函数的零点,我知道如何使用它来计算一个数的平方根,但我看不出这些代码如何执行它.

有人可以向我解释一下吗？

Answer 1

鉴于牛顿迭代 ,在源代码中看到这一点应该很直接.

 __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
 );

确切地说,该算法用于反平方根.

请注意,这仍然无法提供完全准确的结果. rsqrtps使用NR迭代可以得到近23位的精度,而sqrtps最后一位的正确舍入则为24位.

如果要将结果截断为整数,则精度有限是一个问题. (int)4.99999是4.另外,x == 0.0如果使用sqrt(x) ~= x * sqrt(x),请注意这个案例,因为0 * +Inf = NaN.